module user

   !  SUBROUTINES
   !
   !  1) set_wall_temperature
   !  2) set_cp_and_gamma
   !  3) set_variable_cp_and_gamma
   !  4) set_variable_cp_and_gamma_post_proc
   !  5) set_laminar_viscosity_at_nodes
   !  6) set_variable_laminar_viscosity_at_nodes
   !  7) set_variable_laminar_viscosity_at_nodes_post_proc
   !  8) set_thermal_conductivity_at_nodes
   !  9) set_variable_thermal_conductivity_at_nodes
   ! 10) set_variable_thermal_conductivity_at_nodes_post_proc
   ! 11) get_laminar_viscosity_at_faces
   ! 12) get_thermal_conductivity_at_faces
   ! 13) set_bcu
   ! 14) set_bcv
   ! 15) set_bcT
   ! 16) set_bcp
   ! 17) get_uin_vin_pin_Tin_Mw
   ! 18) get_plin_and_p_fictitious
   ! 19) get_u_v_extrapolation_to_fictitious
   ! 20) get_extrapolation_to_corners
   ! 21) get_boundary_simplec_coefficients
   ! 22) get_Uce_Vcn_at_boundary_faces
   ! 23) get_Uce_Vcn_at_boundary_faces_with_pl
   ! 24) get_isentropic_mass_flow
   ! 25) get_mach_area
   ! 26) get_initial_guess
   ! 27) get_boundary_nodes
   ! 28) arbitrary_contour
   ! 29) conical_contour
   ! 30) parabolic_contour
   ! 31) bell_contour
   ! 32) get_spline_coefficients
   ! 33) get_a5d_b_rescaling
   ! 34) get_a9d_b_rescaling

   use coefficients
   use solvers
   use depp_interface

   implicit none

contains

   !============================================================================

   subroutine set_wall_temperature(nx, Tw_cte, Twall) ! Last one is output
      implicit none
      integer, intent(in) :: nx           !< Number of volumes in csi direction (real + fictitious)
      real(8), intent(in) :: Tw_cte       !< constant wall temperature
      real(8), intent(out) :: Twall(nx)   !< Wall temperature

      Twall = Tw_cte

   end subroutine set_wall_temperature

   !============================================================================

   subroutine set_cp_and_gamma(nx, ny, Rg, cp_cte, cp, gcp) ! Last two are output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: Rg  !< Perfect gas constant
      real(8), intent(in) :: cp_cte       !< constant specific heat
      real(8), intent(out) :: cp(nx*ny)   !< Specific heat at const pressure
      real(8), intent(out) :: gcp(nx*ny)  !< gcp = gamma = Cp/Cv at center of CV P

      cp = cp_cte

      gcp = cp/(cp - Rg)

   end subroutine set_cp_and_gamma

   !============================================================================

   subroutine set_variable_cp_and_gamma(folder_input, nx, ny, specie, Rg, T, &
         cp, gcp) ! Last two are output
      implicit none
      character(200), intent(in) :: folder_input   !< input folder name
      integer, intent(in) :: nx     !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: specie !< chemical specie (1 = H20; 2 = AIR)
      real(8), intent(in) :: Rg     !< Perfect gas constant
      real(8), intent(in) :: T(nx*ny)     !< Temperature
      real(8), intent(out) :: cp(nx*ny)   !< Specific heat at const pressure
      real(8), intent(out) :: gcp(nx*ny)  !< gcp = gamma = Cp/Cv at center of CV P

      integer :: np     !< index of the control volume P
      integer :: i, j   !< counters
      integer :: l      !< temperature level (1=T > 1000; 2=T < 1000)
      integer :: m      !< coefficient counter
      integer :: k      !< counters
      integer :: nm     !< number of the measures
      real(8) :: p1     !< weights for the average
      real(8) :: p2     !< weights for the average
      real(8) :: a_cp(5,2)   !< Coefficient matrix of the interpolation function for the specific heat
      real(8), allocatable, dimension(:,:) :: measure

      if (specie == 1) then

         ! One-species, variable properties - Species: H2O

         open(501, file = trim(folder_input) // 'properties01.dat')
         read(501,*)
         read(501,*)
         read(501,*)
         read(501,*)

         do l = 1, 2
            do m = 1, 5
               read(501,*) a_cp(m,l)
            end do
            read(501,*)
            read(501,*)
            read(501,*)
         end do

         close(501)

         do i = 2, nx-1

            do j = 2, ny-1

               np = (j - 1)*nx + i

               if (T(np) >= 1.d3) then
                  l = 1
               else
                  l = 2
               end if

               cp(np) = (a_cp(1,l) + a_cp(2,l)*T(np) + a_cp(3,l)*(T(np)**2)&
                  + a_cp(4,l)*(T(np)**3) + a_cp(5, l)*(T(np)**4))*Rg

            end do

         end do

      else if (specie == 2) then

         ! One-species, variable properties - Species: AIR

         open(501, file = trim(folder_input) // 'properties_air.dat')
         read(501,*)
         read(501,*) nm
         read(501,*)
         read(501,*)

         allocate(measure(nm,8))

         do k = 1, nm
            read(501,*) measure(k, :)
         end do

         close(501)

         do i = 2, nx-1

            do j = 2, ny-1

               np = (j - 1)*nx + i
               l = 1

               do k = 1, nm
                  if (measure(k,1) < T(np)) l = k
               end do

               if (T(np) > measure(1,1) .and. T(np) < measure(nm,1)) then
                  p1 = measure(l,1) - T(np)
                  p2 = T(np) - measure(l+1,1)
                  cp(np) = (p2*measure(l,3) + p1*measure(l+1,3))/(p1 + p2)
               else
                  cp(np) = measure(l,3)
               end if

            end do

         end do

      else

         write(*,*)
         write(*,*) "ERROR: invalid value to SPECIE variable."
         write(*,*)
         stop

      end if

      gcp = cp/(cp - Rg)

   end subroutine set_variable_cp_and_gamma

   !============================================================================

   subroutine set_variable_cp_and_gamma_post_proc(folder_input, nx, ny, &
         specie, Rg, T, cp, gcp) ! Last two are output
      implicit none
      character(200), intent(in) :: folder_input   !< input folder name
      integer, intent(in) :: nx     !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: specie !< chemical specie (1=H20; 2=AIR)
      real(8), intent(in) :: Rg     !< Perfect gas constant
      real(8), intent(in) :: T(nx*ny)     !< Temperature
      real(8), intent(out) :: cp(nx*ny)   !< Specific heat at const pressure
      real(8), intent(out) :: gcp(nx*ny)  !< gcp = gamma = Cp/Cv at center of CV P

      integer :: np     !< Index of the control volume P
      integer :: i, j   !< Counters
      integer :: l      !< Temperature level (1=T > 1000; 2=T < 1000)
      integer :: m      !< Coefficient counter
      integer :: k      !< Counters
      integer :: nm     !< number of the measures
      real(8) :: p1     !< weights for the average
      real(8) :: p2     !< weights for the average
      real(8) :: a_cp(5,2) !< Coefficient matrix of the interpolation function for the specific heat
      real(8), allocatable, dimension(:,:) :: measure

      if (specie == 1) then

         ! One-species, variable properties - Species: H2O

         open(501, file = trim(folder_input) // 'properties01.dat')
         read(501,*)
         read(501,*)
         read(501,*)
         read(501,*)

         do l = 1, 2
            do m = 1, 5
               read(501,*) a_cp(m, l)
            end do
            read(501,*)
            read(501,*)
            read(501,*)
         end do

         close(501)

         do i = 1, nx

            do j = 1, ny

               np = (j - 1)*nx + i

               if (T(np) >= 1.d3) then
                  l = 1
               else
                  l = 2
               end if

               cp(np) = (a_cp(1,l) + a_cp(2,l)*T(np) + a_cp(3,l)*(T(np)**2) &
                  + a_cp(4,l)*(T(np)**3) + a_cp(5,l)*(T(np)**4))*Rg

            end do

         end do

      else if (specie == 2) then

         ! One-species, variable properties - Species: AIR

         open(501, file = trim(folder_input) // 'properties_air.dat')
         read(501,*)
         read(501,*) nm
         read(501,*)
         read(501,*)

         allocate(measure(nm,8))

         do k = 1, nm
            read(501,*) measure(k, :)
         end do

         close(501)

         do i = 1, nx

            do j = 1, ny

               np = (j - 1)*nx + i
               l = 1

               do k = 1, nm
                  if (measure(k,1) < T(np)) l = k
               end do

               if (T(np) > measure(1,1) .and. T(np) < measure(nm,1)) then
                  p1 = measure(l,1) - T(np)
                  p2 = T(np) - measure(l+1,1)
                  cp(np) = (p2*measure(l,3) + p1*measure(l+1,3))/(p1 + p2)
               else
                  cp(np) = measure(l,3)
               end if

            end do

         end do

      else

         write(*,*)
         write(*,*) "ERROR: invalid value to SPECIE variable."
         write(*,*)
         stop

      end if

      gcp = cp/(cp - Rg)

   end subroutine set_variable_cp_and_gamma_post_proc

   !============================================================================

   subroutine set_laminar_viscosity_at_nodes(nx, ny, visc_cte, vlp) ! Last one is output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: visc_cte     !< constant dynamic viscosity
      real(8), intent(out) :: vlp(nx*ny)  !< Laminar viscosity at center of volume P

      vlp = visc_cte

   end subroutine set_laminar_viscosity_at_nodes

   !============================================================================

   subroutine set_variable_laminar_viscosity_at_nodes(folder_input, nx, ny, &
         specie, T, vlp) ! Last one is output
      implicit none
      character(200), intent(in) :: folder_input   !< input folder name
      integer, intent(in) :: nx     !< number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     !< number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: specie !< chemical specie (1=H20; 2=AIR)
      real(8), intent(in) :: T(nx*ny)     !< Temperature
      real(8), intent(out) :: vlp(nx*ny)  !< Laminar viscosity at center of volume P

      integer :: np     !< Index of the control volume P
      integer :: i, j   !< Counters
      integer :: l      !< Temperature level (1=T > 1000; 2=T < 1000)
      integer :: m      !< Coefficient counter
      integer :: k      !< Counters
      integer :: nm     !< number of the measures
      real(8) :: p1     !< weights for the average
      real(8) :: p2     !< weights for the average
      real(8) :: a_vlp(4,2)   !< Coefficient matrix of the interpolation function for the conductivity
      real(8), allocatable, dimension(:,:) :: measure

      if (specie == 1) then

         ! One-species, variable properties

         open(501, file = trim(folder_input) // 'properties02.dat')
         read(501,*)
         read(501,*)

         do l = 1, 2
            read(501,*) (a_vlp(m,l), m = 1,4)
         end do

         close(501)

         do i = 2, nx-1

            do j = 2, ny-1

               np = (j - 1)*nx + i

               if (T(np) >= 1.d3) then
                  l = 1
               else
                  l = 2
               end if

               vlp(np) = dexp((a_vlp(1,l)*dlog(T(np)) + a_vlp(2,l)/T(np) &
                  + a_vlp(3,l)/(T(np)**2) + a_vlp(4,l)))*1.d-7

            end do

         end do

      else if (specie == 2) then

         ! One-species, variable properties - Species: AIR

         open(501, file = trim(folder_input) // 'properties_air.dat')
         read(501,*)
         read(501,*) nm
         read(501,*)
         read(501,*)

         allocate(measure(nm,8))

         do k = 1, nm
            read(501,*) measure(k, :)
         end do

         close(501)

         do i = 2, nx-1

            do j = 2, ny-1

               np = (j - 1)*nx + i
               l = 1

               do k = 1, nm
                  if (measure(k,1) < T(np)) l = k
               end do

               if (T(np) > measure(1,1) .and. T(np) < measure(nm,1)) then
                  p1 = measure(l,1) - T(np)
                  p2 = T(np) - measure(l+1,1)
                  vlp(np) = (p2*measure(l,4) + p1*measure(l+1,4))/(p1 + p2)
               else
                  vlp(np) = measure(l,4)
               end if

            end do

         end do

      else

         write(*,*)
         write(*,*) "ERROR: invalid value to SPECIE variable."
         write(*,*)
         stop

      end if

   end subroutine set_variable_laminar_viscosity_at_nodes

   !============================================================================

   subroutine set_variable_laminar_viscosity_at_nodes_post_proc(folder_input, &
         nx, ny, specie, T, vlp) ! Last variable is output
      implicit none
      character(200), intent(in) :: folder_input   !< input folder name
      integer, intent(in) :: nx           !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny           !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: specie       !< chemical specie (1 = H20; 2 = AIR)
      real(8), intent(in) :: T(nx*ny)     !< Temperature
      real(8), intent(out) :: vlp(nx*ny)  !< Laminar viscosity at center of volume P

      integer :: np     !< Index of the control volume P
      integer :: i, j   !< Counters
      integer :: l      !< Temperature level (1=T > 1000; 2=T < 1000)
      integer :: m      !< Coefficient counter
      integer :: k      !< counters
      integer :: nm     !< number of the measures
      real(8) :: p1     !< weights for the average
      real(8) :: p2     !< weights for the average
      real(8) :: a_vlp(4,2)   !< Coefficient matrix of the interpolation function for the conductivity
      real(8), allocatable, dimension(:,:) :: measure

      if (specie == 1) then

         ! One-species, variable properties

         open(501, file = trim(folder_input) // 'properties02.dat')
         read(501,*)
         read(501,*)

         do l = 1, 2
            read(501,*) (a_vlp(m, l), m = 1,4)
         end do

         close(501)

         do i = 1, nx

            do j = 1, ny

               np = (j - 1)*nx + i

               if (T(np) >= 1.d3) then
                  l = 1
               else
                  l = 2
               end if

               vlp(np) = dexp((a_vlp(1,l)*dlog(T(np)) + a_vlp(2,l)/T(np) &
                  + a_vlp(3,l)/(T(np)**2) + a_vlp(4,l)))*1.d-7

            end do

         end do

      else if (specie == 2) then

         ! One-species, variable properties - Species: AIR

         open(501, file = trim(folder_input) // 'properties_air.dat')
         read(501,*)
         read(501,*) nm
         read(501,*)
         read(501,*)

         allocate(measure(nm,8))

         do k = 1, nm
            read(501,*) measure(k, :)
         end do

         close(501)

         do i = 1, nx

            do j = 1, ny

               np = (j - 1)*nx + i
               l = 1

               do k = 1, nm
                  if (measure(k,1) < T(np)) l = k
               end do

               if (T(np) > measure(1,1) .and. T(np) < measure(nm,1)) then
                  p1 = measure(l,1) - T(np)
                  p2 = T(np) - measure(l+1,1)
                  vlp(np) = (p2*measure(l,4) + p1*measure(l+1,4))/(p1 + p2)
               else
                  vlp(np) = measure(l,4)
               end if

            end do

         end do

      else

         write(*,*)
         write(*,*) "ERROR: invalid value to SPECIE variable."
         write(*,*)
         stop

      end if

   end subroutine set_variable_laminar_viscosity_at_nodes_post_proc

   !============================================================================

   subroutine set_thermal_conductivity_at_nodes(nx, ny, k_cte, kp) ! Last one is output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: k_cte        !< Constant thermal conductivity
      real(8), intent(out) :: kp(nx*ny)   !< Thermal conductivity at center of volume P

      kp = k_cte

   end subroutine set_thermal_conductivity_at_nodes

   !============================================================================

   subroutine set_variable_thermal_conductivity_at_nodes(folder_input, nx, ny, &
         specie, T, kp) ! Last one is output
      implicit none
      character(200), intent(in) :: folder_input   !< input folder name
      integer, intent(in) :: nx     !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: specie !< chemical specie (1 = H20; 2 = AIR)
      real(8), intent(in) :: T(nx*ny)     !< Temperature
      real(8), intent(out) :: kp(nx*ny)   !< Thermal conductivity at center of volume P

      integer :: np     !< Index of the control volume P
      integer :: i, j   !< Counters
      integer :: l      !< Temperature level (1=T > 1000; 2=T < 1000)
      integer :: m      !< Coefficient counter
      integer :: k      !< Counters
      integer :: nm     !< Number of the measures
      real(8) :: p1     !< Weights for the average
      real(8) :: p2     !< Weights for the average
      real(8) :: a_kp(4,2) !< Coefficient matrix of the interpolation function for the conductivity
      real(8), allocatable, dimension(:,:) :: measure

      if (specie == 1) then

         ! One-species, variable properties - Species: H2O

         open (501, file = trim(folder_input) // 'properties02.dat')
         read(501,*)
         read(501,*)
         read(501,*)
         read(501,*)

         do l = 1, 2
            read(501,*) (a_kp(m, l), m = 1,4)
         end do

         close(501)

         do i = 2, nx-1

            do j = 2, ny-1

               np = (j - 1)*nx + i

               if (T(np) >= 1.d3) then
                  l = 1
               else
                  l = 2
               end if

               kp(np) = dexp((a_kp(1,l)*dlog(T(np)) + a_kp(2,l)/T(np) &
                  + a_kp(3,l)/T(np)**2 + a_kp(4,l)))*1.d-4

            end do

         end do

      else if (specie == 2) then

         ! One-species, variable properties - Species: AIR

         open(501, file = trim(folder_input) // 'properties_air.dat')
         read(501,*)
         read(501,*) nm
         read(501,*)
         read(501,*)

         allocate(measure(nm,8))

         do k = 1, nm
            read(501,*) measure(k, :)
         end do

         close(501)

         do i = 2, nx-1

            do j = 2, ny-1

               np = (j - 1)*nx + i
               l = 1

               do k = 1, nm
                  if (measure(k,1) < T(np)) l = k
               end do

               if (T(np) > measure(1,1) .and. T(np) < measure(nm,1)) then
                  p1 = measure(l,1) - T(np)
                  p2 = T(np) - measure(l + 1,1)
                  kp(np) = (p2*measure(l,6) + p1*measure(l + 1,6))/(p1 + p2)
               else
                  kp(np) = measure(l,6)
               end if

            end do

         end do

      else

         write(*,*)
         write(*,*) "ERROR: invalid value to SPECIE variable."
         write(*,*)
         stop

      end if

   end subroutine set_variable_thermal_conductivity_at_nodes

   !============================================================================

   subroutine set_variable_thermal_conductivity_at_nodes_post_proc( &
         folder_input, nx, ny, specie, T, kp) ! Last one is output
      implicit none
      character(200), intent(in) :: folder_input   !< input folder name
      integer, intent(in) :: nx     !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: specie !< Chemical specie (1 = H20; 2 = AIR)
      real(8), intent(in) :: T(nx*ny)     !< Temperature
      real(8), intent(out) :: kp(nx*ny)   !< Thermal conductivity at center of volume P

      integer :: np     !< Index of the control volume P
      integer :: i, j   !< Counters
      integer :: l      !< Temperature level (1=T > 1000; 2=T < 1000)
      integer :: m      !< Coefficient counter
      integer :: k      !< Counters
      integer :: nm     !< Number of the measures
      real(8) :: p1     !< Weights for the average
      real(8) :: p2     !< Weights for the average
      real(8) :: a_kp(4,2) !< Coefficient matrix of the interpolation function for the conductivity
      real(8), allocatable, dimension(:,:) :: measure

      if (specie == 1) then

         ! One-species, variable properties - Species: H2O

         open (501, file = trim(folder_input) // 'properties02.dat')
         read(501,*)
         read(501,*)
         read(501,*)
         read(501,*)

         do l = 1, 2
            read(501,*) (a_kp(m, l), m = 1,4)
         end do

         close(501)

         do i = 1, nx

            do j = 1, ny

               np = (j - 1)*nx + i

               if (T(np) >= 1.d3) then
                  l = 1
               else
                  l = 2
               end if

               kp(np) = dexp((a_kp(1,l)*dlog(T(np)) + a_kp(2,l)/T(np) &
                  + a_kp(3,l)/T(np)**2 + a_kp(4,l)))*1.d-4

            end do

         end do

      else if (specie == 2) then

         ! One-species, variable properties - Species: AIR

         open(501, file = trim(folder_input) // 'properties_air.dat')
         read(501,*)
         read(501,*) nm
         read(501,*)
         read(501,*)

         allocate(measure(nm,8))

         do k = 1, nm
            read(501,*) measure(k, :)
         end do

         close(501)

         do i = 1, nx

            do j = 1, ny

               np = (j - 1)*nx + i
               l = 1

               do k = 1, nm
                  if (measure(k,1) < T(np)) l = k
               end do

               if (T(np) > measure(1,1) .and. T(np) < measure(nm,1)) then
                  p1 = measure(l,1) - T(np)
                  p2 = T(np) - measure(l+1,1)
                  kp(np) = (p2*measure(l,6) + p1*measure(l+1,6))/(p1 + p2)
               else
                  kp(np) = measure(l,6)
               end if

            end do

         end do

      else

         write(*,*)
         write(*,*) "ERROR: invalid value to SPECIE variable."
         write(*,*)
         stop

      end if

   end subroutine set_variable_thermal_conductivity_at_nodes_post_proc

   !============================================================================

   ! CAUTION: In order to apply correctly this function, laminar viscosity vlp must be known
   ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
   subroutine get_laminar_viscosity_at_faces(nx, ny, vlp, vle, vln) ! Last two are output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: vlp(nx*ny)   !< Laminar viscosity at center of volume P
      real(8), intent(out) :: vle(nx*ny)  !< Laminar viscosity at center of face east
      real(8), intent(out) :: vln(nx*ny)  !< Laminar viscosity at center of face north

      integer :: i, j, np, npe, npn

      do j = 2, ny-1

         do i = 1, nx-1

            np = nx*(j - 1) + i
            npe = np + 1

            vle(np) = 2.d0*vlp(np)*vlp(npe)/(vlp(np) + vlp(npe))

         end do

      end do

      do i = 2, nx-1

         do j = 1, ny-1

            np = nx*(j - 1) + i
            npn = np + nx

            vln(np) = 2.d0*vlp(np)*vlp(npn)/(vlp(np) + vlp(npn))

         end do

      end do

   end subroutine get_laminar_viscosity_at_faces

   !============================================================================

   ! CAUTION: In order to apply correctly this function, thermal conductivity kp must be known
   ! in real and fictitious volumes (except CV of corners SW, SE, NW and NE)
   subroutine get_thermal_conductivity_at_faces(nx, ny, kp, ke, kn) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: kp(nx*ny)    !< Thermal conductivity at center of volume P
      real(8), intent(out) :: ke(nx*ny)   !< Thermal conductivity at center of face east
      real(8), intent(out) :: kn(nx*ny)   !< Thermal conductivity at center of face north

      integer :: i, j, np, npe, npn

      do j = 2, ny-1

         do i = 1, nx-1

            np = nx*(j - 1) + i
            npe = np + 1

            ke(np) = 2.d0*kp(np)*kp(npe)/(kp(np) + kp(npe))

         end do

      end do

      do i = 2, nx-1

         do j = 1, ny-1

            np = nx*(j - 1) + i
            npn = np + nx

            kn(np) = 2.d0*kp(np)*kp(npn)/(kp(np) + kp(npn))

         end do

      end do

   end subroutine get_thermal_conductivity_at_faces

   !============================================================================

   subroutine set_bcu(nx, ny, modvis, u, betan, gamman, au, bu) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny              !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: modvis          !< Viscosity model (0=Euler; 1=Navier-Stokes)
      real(8), intent(in) :: u(nx*ny)        !< Cartesian velocity of the last iteraction
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: au(nx*ny,9)  !< Coefficients of the linear system for u
      real(8), intent(inout) :: bu(nx*ny)    !< Source vector of the linear system

      integer :: i, j, np
      real(8) :: ubn(nx)

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         au(np,:) = 0.d0
         au(np,6) = -1.d0
         au(np,5) = 1.d0
         bu(np) = 0.d0
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         au(np,:) = 0.d0
         au(np,4) = -1.d0
         au(np,5) = 1.d0
         bu(np) = 0.d0
      end do

      ! South boundary
      call set_null_normal_grad_south(nx, ny, gamman, betan, au, bu)

      ! North boundary
      if (modvis == 1) then ! Navier Stokes

         ubn = 0.d0
         call set_prescribed_condition_north(nx, ny, ubn, au, bu)

      else ! Euler

         call set_null_normal_grad_north(nx, ny, gamman, betan, au, bu)

      end if

      ! Corners (extrapolation)
      call set_coef_extrapolation_to_corners(nx, ny, u, au, bu)

   end subroutine set_bcu

   !============================================================================

   subroutine set_bcv(nx, ny, modvis, v, vbw, betan, gamman, av, bv) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny              !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: modvis          !< Viscosity model (0=Euler; 1=Navier-Stokes)
      real(8), intent(in) :: v(nx*ny)        !< Cartesian velocity of the last iteraction
      real(8), intent(in) :: vbw(ny)         !< v over the faces of the west boundary (m/s)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(out) :: av(nx*ny,9)    !< Coefficients of the linear system for v
      real(8), intent(out) :: bv(nx*ny)      !< Source vector of the linear system

      integer :: i, j, np
      real(8) :: vbn(nx)   !< v over the faces of the north boundary (m/s)
      real(8) :: vbs(nx)   !< v over the faces of the south boundary (m/s)

      ! West boundary
      call set_prescribed_condition_west(nx, ny, vbw, av, bv)

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         av(np,:) = 0.d0
         av(np,4) = -1.d0
         av(np,5) = 1.d0
         bv(np) = 0.d0
      end do

      ! South boundary
      vbs = 0.d0
      call set_prescribed_condition_south(nx, ny, vbs, av, bv)

      ! North boundary
      if (modvis == 1) then ! Navier-Stokes

         vbn = 0.d0
         call set_prescribed_condition_north(nx, ny, vbn, av, bv)

      else ! Euler

         call set_null_normal_grad_north(nx, ny, gamman, betan, av, bv)

      end if

      ! Corners (extrapolation)
      call set_coef_extrapolation_to_corners(nx, ny, v, av, bv)

   end subroutine set_bcv

   !============================================================================

   ! CAUTION: Twall must be difined in the fictitious volumes too
   subroutine set_bcT(nx, ny, ccTw, Tin, Twall, T, betan, gamman, at, bt) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny              !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: ccTw            !< Boundary condition (0=adiabatic; 1=Dirichlet)
      real(8), intent(in) :: Tin(ny)         !< Temperature in the entrance
      real(8), intent(in) :: Twall(nx)       !< Wall temperature
      real(8), intent(in) :: T(nx*ny)        !< Temperature of the last iteraction
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(out) :: aT(nx*ny,9)    !< Coefficients of the linear system
      real(8), intent(out) :: bt(nx*ny)      !< Source vector of the linear system

      integer :: i, j, np

      ! West boundary
      call set_prescribed_condition_west(nx, ny, Tin, at, bt)

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         aT(np,:) = 0.d0
         aT(np,4) = -1.d0
         aT(np,5) = 1.d0
         bT(np) = 0.d0
      end do

      ! South boundary
      call set_null_normal_grad_south(nx, ny, gamman, betan, at, bt)

      ! North boundary
      if (ccTw == 0) then ! Adiabatic boundary condition

         call set_null_normal_grad_north(nx, ny, gamman, betan, at, bt)

      else ! Temperature prescribed boundary condition

         call set_prescribed_condition_north(nx, ny, Twall, at, bt)

      end if

      ! Corners (extrapolation)
      call set_coef_extrapolation_to_corners(nx, ny, T, at, bt)

   end subroutine set_bcT

   !============================================================================

   !> Calculates the coefficients and source of the linear system of the pressure deviation
   !! for the fictitious volumes based on the boundary conditions
   subroutine set_bcp(nx, ny, ap, bp) ! Output: last two
      implicit none
      integer, intent(in) :: nx           !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny           !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(out) :: ap(nx*ny,5) !< Coefficients of the linear system for pl (m.s)
      real(8), intent(out) :: bp(nx*ny)   !< Source of the linear system for pl (kg/s)

      ! Auxliary variables
      integer :: i, j, np, iaux

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         ap(np,:) = 0.d0
         ap(np,3) = 1.d0
         ap(np,4) = 1.d0
         bp(np) = 0.d0
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         ap(np,:) = 0.d0
         ap(np,2) = -1.d0
         ap(np,3) = 1.d0
         bp(np) = 0.d0
      end do

      ! South boundary
      j = 1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         ap(np,:) = 0.d0
         ap(np,3) = 1.d0
         ap(np,5) = -1.d0
         bp(np) = 0.d0
      end do

      ! North boundary
      j = ny
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         ap(np,:) = 0.d0
         ap(np,1) = -1.d0
         ap(np,3) = 1.d0
         bp(np) = 0.d0
      end do

      ! Corners

      ! SW
      i = 1
      j = 1
      np = nx*(j - 1) + i
      ap(np,:) = 0.d0
      ap(np,3) = 1.d0
      bp(np) = 0.d0

      !SE
      i = nx
      j = 1
      np = nx*(j - 1) + i
      ap(np,:) = 0.d0
      ap(np,3) = 1.d0
      bp(np) = 0.d0

      ! NW
      i = 1
      j = ny
      np = nx*(j - 1) + i
      ap(np,:) = 0.d0
      ap(np,3) = 1.d0
      bp(np) = 0.d0

      ! NE
      i = nx
      j = ny
      np = nx*(j - 1) + i
      ap(np,:) = 0.d0
      ap(np,3) = 1.d0
      bp(np) = 0.d0

   end subroutine set_bcp

   !============================================================================

   !> Extrapolates the field u of the real volumes to the fictitious volumes
   !! based on the boundary conditions.
   subroutine get_u_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, u) ! Last is inout
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      integer, intent(in) :: modvis          !< Viscosity model (0=Euler; 1=Navier-Stokes)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(inout) :: u(nx*ny)     !< Temperature of the last iteraction

      ! Auxiliary variables
      integer :: i, j, iaux, np, npe, nps, npw

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         u(np) = u(npe)
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         u(np) = u(npw)
      end do

      ! South boundary
      call get_extrapolation_null_normal_grad_south(nx, ny, gamman, betan, u)

      ! North boundary
      if (modvis == 0) then ! Euler

         call get_extrapolation_null_normal_grad_north(nx, ny, gamman, betan, u)

      else ! Navier-Stokes

         j = ny
         iaux = nx*(j - 1)
         do i = 2, nx-1
            np = iaux + i
            nps = np - nx
            u(np) = -u(nps)
         end do

      end if

      call get_extrapolation_to_corners(nx, ny, u) ! InOutput: last one

   end subroutine get_u_extrapolation_to_fictitious

   !============================================================================

   !> Extrapolates the field v of the real volumes to the fictitious volumes
   !! based on the boundary conditions.
   subroutine get_v_extrapolation_to_fictitious(nx, ny, modvis, betan, gamman, v) ! Last is inout
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      integer, intent(in) :: modvis          !< Viscosity model (0=Euler; 1=Navier-Stokes)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(inout) :: v(nx*ny)     !< Temperature of the last iteraction

      ! Auxiliary variables
      integer :: i, j, iaux, np, npe, npn, nps, npw

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         v(np) = -v(npe)
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         v(np) = v(npw)
      end do

      ! South boundary
      j = 1
      do i = 2, nx-1
         np = nx*(j - 1) + i
         npn = np + nx
         v(np) = -v(npn)
      end do

      ! North boundary
      if (modvis == 0) then ! Euler

         call get_extrapolation_null_normal_grad_north(nx, ny, gamman, betan, v)

      else ! Navier-Stokes

         j = ny
         iaux = nx*(j - 1)
         do i = 2, nx-1
            np = iaux + i
            nps = np - nx
            v(np) = -v(nps)
         end do

      end if

      call get_extrapolation_to_corners(nx, ny, v) ! InOutput: last one

   end subroutine get_v_extrapolation_to_fictitious

   !============================================================================

   !> Extrapolates the field p of the real volumes to the fictitious volumes
   !! based on the boundary conditions.
   subroutine get_p_extrapolation_to_fictitious(nx, ny, pbw, p) ! Last is inout
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: pbw(ny)         !< p over the faces of the west boundary (Pa)
      real(8), intent(inout) :: p(nx*ny)     !< F at center of volume P (K)

      ! Auxiliary variables
      integer :: i, j, np, npe, npw, nps, npn

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         p(np) = 2.d0*pbw(j) - p(npe)
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         p(np) = p(npw)
      end do

      ! South boundary
      j = 1
      do i = 2, nx-1
         np = nx*(j - 1) + i
         npn = np + nx
         p(np) = p(npn)
      end do

      ! North boundary
      j = ny
      do i = 2, nx-1
         np = nx*(j - 1) + i
         nps = np - nx
         p(np) = p(nps)
      end do

      call get_extrapolation_to_corners(nx, ny, p) ! InOutput: last one

   end subroutine get_p_extrapolation_to_fictitious

   !============================================================================

   !> \brief Calculates the coefficients for linear system assuming that
   !! the variables are constant along the streamline
   subroutine set_streamlined_exit_east(nx, ny, Ucbe, Vcbe, af, bf)
      implicit none
      integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: Ucbe(ny)        !< Uc over the faces of the east boundary (m2/s)
      real(8), intent(in) :: Vcbe(ny)        !< Vc over the faces of the east boundary (m2/s)
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables

      integer :: i, j, np

      i = nx

      ! Calculating coefficients for the real volume of the corner SE

      j = 2

      np = nx*(j - 1) + i

      af(np,:) = 0.d0

      af(np,4) = -1.d0 - 0.5d0*Vcbe(j)/Ucbe(j)  ! West
      af(np,5) = 1.d0 - 0.5d0*Vcbe(j)/Ucbe(j)   ! P
      af(np,7) = 0.5d0*Vcbe(j)/Ucbe(j)          ! North-west
      af(np,8) = 0.5d0*Vcbe(j)/Ucbe(j)          ! North

      bf(np) = 0.d0

      ! Calculating coefficients for the real volumes of the east boundary, except SE and NE

      do j = 3, ny-2

         np = nx*(j - 1) + i

         af(np,:) = 0.d0

         af(np,1) = -0.25d0*Vcbe(j)/Ucbe(j)  ! South-west
         af(np,2) = -0.25d0*Vcbe(j)/Ucbe(j)  ! South
         af(np,4) = -1.d0                    ! West
         af(np,5) = 1.d0                     ! P
         af(np,7) = 0.25d0*Vcbe(j)/Ucbe(j)   ! North-west
         af(np,8) = 0.25d0*Vcbe(j)/Ucbe(j)   ! North

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the real volume of the corner NE

      j = ny - 1

      np = nx*(j - 1) + i

      af(np,:) = 0.d0

      af(np,1) = -0.5d0*Vcbe(j)/Ucbe(j)         ! South-west
      af(np,2) = -0.5d0*Vcbe(j)/Ucbe(j)         ! South
      af(np,4) = -1.d0 + 0.5d0*Vcbe(j)/Ucbe(j)  ! West
      af(np,5) = 1.d0 + 0.5d0*Vcbe(j)/Ucbe(j)   ! P

      bf(np) = 0.d0

   end subroutine set_streamlined_exit_east

   !============================================================================

   !> \brief Calculates the coefficients for linear system assuming that
   !! the variables are constant along the streamline
   subroutine set_streamlined_exit_east_5d(nx, ny, Ucbe, Vcbe, af, bf)
      implicit none
      integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: Ucbe(ny)        !< Uc over the faces of the east boundary (m2/s)
      real(8), intent(in) :: Vcbe(ny)        !< Vc over the faces of the east boundary (m2/s)
      real(8), intent(inout) :: af(nx*ny,5)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables

      integer :: i, j, np

      i = nx

      ! Calculating coefficients for the real volume of the corner SE

      j = 2

      np = nx*(j - 1) + i

      af(np,:) = 0.d0

      af(np,2) = -1.d0                    ! West
      af(np,3) = 1.d0 - Vcbe(j)/Ucbe(j)   ! P
      af(np,5) = Vcbe(j)/Ucbe(j)          ! North

      bf(np) = 0.d0

      ! Calculating coefficients for the real volumes of the east boundary, except SE and NE

      do j = 3, ny-2

         np = nx*(j - 1) + i

         af(np,:) = 0.d0

         af(np,1) = -0.5d0*Vcbe(j)/Ucbe(j)   ! South
         af(np,2) = -1.d0                    ! West
         af(np,3) = 1.d0                     ! P
         af(np,5) = 0.5d0*Vcbe(j)/Ucbe(j)    ! North

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the real volume of the corner NE

      j = ny - 1

      np = nx*(j - 1) + i

      af(np,:) = 0.d0

      af(np,1) = -Vcbe(j)/Ucbe(j)         ! South
      af(np,2) = -1.d0                    ! West
      af(np,3) = 1.d0 + Vcbe(j)/Ucbe(j)   ! P

      bf(np) = 0.d0

   end subroutine set_streamlined_exit_east_5d

   !============================================================================

   !> \brief Calculates the coefficients for linear system assuming that
   !! the variables are constant along the streamline
   subroutine get_extrapolation_streamlined_east(nx, ny, Ucbe, Vcbe, f)
      implicit none
      integer, intent(in) :: nx   !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny   !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: Ucbe(ny)        !< Uc over the faces of the east boundary (m2/s)
      real(8), intent(in) :: Vcbe(ny)        !< Vc over the faces of the east boundary (m2/s)
      real(8), intent(inout) :: F(nx*ny)     !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, np, npw, npsw, npnw
      real(8) :: faux(ny-2)
      real(8) :: bf(ny-2)
      real(8) :: af(ny-2,3)

      i = nx

      ! Calculating coefficients for the real volume of the corner SE

      j = 2

      np = nx*(j - 1) + i
      npw = np - 1
      npnw = npw + nx

      af(j-1,1) = 0.d0                          ! South
      af(j-1,2) = 1.d0 - 0.5d0*Vcbe(j)/Ucbe(j)  ! P
      af(j-1,3) = 0.5d0*Vcbe(j)/Ucbe(j)         ! North

      bf(j-1) = F(npw) - 0.5d0*Vcbe(j)/Ucbe(j)*(F(npnw) - F(npw))

      ! Calculating coefficients for the real volumes of the east boundary, except SE and NE

      do j = 3, ny-2

         np = nx*(j - 1) + i
         npw = np - 1
         npsw = npw - nx
         npnw = npw + nx

         af(j-1,1) = -0.25d0*Vcbe(j)/Ucbe(j) ! South
         af(j-1,2) = 1.d0                    ! P
         af(j-1,3) = 0.25d0*Vcbe(j)/Ucbe(j)  ! North

         bf(j-1) = F(npw) - 0.25d0*Vcbe(j)/Ucbe(j)*(F(npnw) - F(npsw))

      end do

      ! Calculating coefficients for the real volume of the corner NE

      j = ny - 1

      np = nx*(j - 1) + i
      npw = np - 1
      npsw = npw - nx

      af(j-1,1) = -0.5d0*Vcbe(j)/Ucbe(j)         ! South
      af(j-1,2) = 1.d0 + 0.5d0*Vcbe(j)/Ucbe(j)   ! P

      bf(j-1) = F(npw) - 0.5d0*Vcbe(j)/Ucbe(j)*(F(npw) - F(npsw))

      call tdma3d(ny-2, af, bf, faux)

      do j = 2, ny-1

         np = nx*(j - 1) + i

         F(np) = faux(j-1)

      end do

   end subroutine get_extrapolation_streamlined_east

   !============================================================================

   !> Calculates the field over the south boundary based on the null normal gradient boundary condition
   subroutine get_extrapolation_null_normal_grad_north(nx, ny, gamman, betan, f)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: f(nx*ny)     !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np, nps, npsw, npse
      real(8) :: faux(nx-2)
      real(8) :: bf(nx-2)
      real(8) :: af(nx-2,3)

      af = 0.d0

      bf = 0.d0

      j = ny

      iaux = nx*(j - 1)

      ! Calculating coefficients for the fictitious volume of the corner SW

      i = 2

      np = iaux + i
      nps = np - nx
      npse = nps + 1

      af(i-1,1) = 0.d0                                ! West
      af(i-1,2) = 1.d0 + 0.5d0*betan(nps)/gamman(nps) ! P
      af(i-1,3) = -0.5d0*betan(nps)/gamman(nps)       ! East

      bf(i-1) = f(nps) + 0.5d0*betan(nps)*(f(npse) - f(nps))/gamman(nps)

      ! Calculating coefficients for the fictitious volumes of the south boundary, except SW and SE

      do i = 3, nx-2

         np = iaux + i
         nps = np - nx
         npsw = nps - 1
         npse = nps + 1

         af(i-1,1) = 0.25d0*betan(nps)/gamman(nps)    ! West
         af(i-1,2) = 1.d0                             ! P
         af(i-1,3) = -0.25d0*betan(nps)/gamman(nps)   ! East

         bf(i-1) = f(nps) + 0.25d0*betan(nps)*(f(npse) - f(npsw))/gamman(nps)

      end do

      ! Calculating coefficients for the fictitious volume of the corner SE

      i = nx - 1

      np = iaux + i
      nps = np - nx
      npsw = nps - 1

      af(i-1,1) = 0.5d0*betan(nps)/gamman(nps)        ! West
      af(i-1,2) = 1.d0 - 0.5d0*betan(nps)/gamman(nps) ! P
      af(i-1,3) = 0.d0                                ! North

      bf(i-1) = f(nps) + 0.5d0*betan(nps)*(f(nps) - f(npsw))/gamman(nps)

      call tdma3d(nx-2, af, bf, faux)

      do i = 2, nx-1

         np = iaux + i

         F(np) = faux(i-1)

      end do

   end subroutine get_extrapolation_null_normal_grad_north

   !============================================================================

   !> Calculates the field over the south boundary based on the null normal gradient boundary condition
   subroutine get_extrapolation_null_normal_grad_south(nx, ny, gamman, betan, f)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: f(nx*ny)     !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np, npn, npnw, npne
      real(8) :: faux(nx-2)
      real(8) :: bf(nx-2)
      real(8) :: af(nx-2,3)

      j = 1

      iaux = nx*(j - 1)

      ! Calculating coefficients for the fictitious volume of the corner SW

      i = 2

      np = iaux + i
      npn = np + nx
      npne = npn + 1

      af(i-1,1) = 0.d0                                ! West
      af(i-1,2) = -1.d0 + 0.5d0*betan(np)/gamman(np)  ! P
      af(i-1,3) = -0.5d0*betan(np)/gamman(np)         ! East

      bf(i-1) = -f(npn) + 0.5d0*betan(np)*(f(npne) - f(npn))/gamman(np)

      ! Calculating coefficients for the fictitious volumes of the south boundary, except SW and SE

      do i = 3, nx-2

         np = iaux + i
         npn = np + nx
         npnw = npn - 1
         npne = npn + 1

         af(i-1,1) = 0.25d0*betan(np)/gamman(np)   ! West
         af(i-1,2) = -1.d0                         ! P
         af(i-1,3) = -0.25d0*betan(np)/gamman(np)  ! East

         bf(i-1) = -f(npn) + 0.25d0*betan(np)*(f(npne) - f(npnw))/gamman(np)

      end do

      ! Calculating coefficients for the fictitious volume of the corner SE

      i = nx - 1

      np = iaux + i
      npn = np + nx
      npnw = npn - 1

      af(i-1,1) = 0.5d0*betan(np)/gamman(np)          ! West
      af(i-1,2) = -1.d0 - 0.5d0*betan(np)/gamman(np)  ! P
      af(i-1,3) = 0.d0                                ! North

      bf(i-1) = -f(npn) + 0.5d0*betan(np)*(f(npn) - f(npnw))/gamman(np)

      call tdma3d(nx-2, af, bf, faux)

      do i = 2, nx-1

         np = iaux + i

         F(np) = faux(i-1)

      end do

   end subroutine get_extrapolation_null_normal_grad_south

   !============================================================================

   !> Calculates the field over the north boundary based on the prescribed boundary condition
   subroutine set_prescribed_condition_north(nx, ny, fbn, af, bf)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: fbn(nx)         !< Values of the field over the north boundary
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables

      integer :: i, j, iaux, np

      ! Calculating coefficients for the fictitious volumes of the north boundary

      j = ny

      iaux = nx*(j - 1)

      do i = 2, nx-1

         np = iaux + i

         af(np,:) = 0.d0

         af(np,2) = 1.d0

         af(np,5) = 1.d0

         bf(np) = 2.d0*fbn(i)

      end do

   end subroutine set_prescribed_condition_north

   !============================================================================

   !> Calculates the field over the north boundary based on the prescribed boundary condition
   subroutine set_prescribed_condition_south(nx, ny, fbs, af, bf)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: fbs(nx)         !< Values of the field over the north boundary
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np

      ! Calculating coefficients for the fictitious volumes of the north boundary

      j = 1

      iaux = nx*(j - 1)

      do i = 2, nx-1

         np = iaux + i

         af(np,:) = 0.d0

         af(np,2) = 1.d0

         af(np,5) = 1.d0

         bf(np) = 2.d0*fbs(i)

      end do

   end subroutine set_prescribed_condition_south

   !============================================================================

   !> Calculates the field over the west boundary based on the prescribed value boundary condition
   subroutine set_prescribed_condition_west(nx, ny, fbw, af, bf)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: fbw(ny)         !< Values of the field over the west boundary
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables

      integer :: i, j, np

      ! Calculating coefficients for the fictitious volumes of the west boundary

      i = 1

      do j = 2, ny-1

         np = nx*(j - 1) + i

         af(np,:) = 0.d0

         af(np,6) = 1.d0

         af(np,5) = 1.d0

         bf(np) = 2.d0*fbw(j)

      end do

   end subroutine set_prescribed_condition_west

   !============================================================================

   !> Calculates the field over the west boundary based on the prescribed value boundary condition
   subroutine set_prescribed_condition_west_5d(nx, ny, fbw, af, bf)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: fbw(ny)         !< Values of the field over the west boundary
      real(8), intent(inout) :: af(nx*ny,5)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables

      integer :: i, j, np

      ! Calculating coefficients for the fictitious volumes of the west boundary

      i = 1

      do j = 2, ny-1

         np = nx*(j - 1) + i

         af(np,:) = 0.d0

         af(np,4) = 1.d0

         af(np,3) = 1.d0

         bf(np) = 2.d0*fbw(j)

      end do

   end subroutine set_prescribed_condition_west_5d

   !============================================================================

   !> Calculates the field over the north boundary based on the null normal gradient boundary condition
   subroutine set_null_normal_grad_north(nx, ny, gamman, betan, af, bf)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np, nps

      j = ny

      iaux = nx*(j - 1)

      ! Calculating coefficients for the fictitious volume of the corner NW

      i = 2

      np = iaux + i
      nps = np - nx

      af(np,:) = 0.d0

      af(np,2) = -1.d0 + 0.5d0*betan(nps)/gamman(nps) ! South
      af(np,3) = -0.5d0*betan(nps)/gamman(nps)        ! Southeast
      af(np,5) = 1.d0 + 0.5d0*betan(nps)/gamman(nps)  ! P
      af(np,6) = -0.5d0*betan(nps)/gamman(nps)        ! East

      bf(np) = 0.d0

      ! Calculating coefficients for the fictitious volumes of the north boundary, except NW and NE

      do i = 3, nx-2

         np = iaux + i
         nps = np - nx

         af(np,:) = 0.d0

         af(np,1) = 0.25d0*betan(nps)/gamman(nps)  ! South-west
         af(np,2) = -1.d0                          ! South
         af(np,3) = -0.25d0*betan(nps)/gamman(nps) ! Southeast
         af(np,4) = 0.25d0*betan(nps)/gamman(nps)  ! West
         af(np,5) = 1.d0                           ! P
         af(np,6) = -0.25d0*betan(nps)/gamman(nps) ! East

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the fictitious volume of the corner NE

      i = nx - 1

      np = iaux + i
      nps = np - nx

      af(np,:) = 0.d0

      af(np,1) = 0.5d0*betan(nps)/gamman(nps)         ! South-west
      af(np,2) = -1.d0 - 0.5d0*betan(nps)/gamman(nps) ! South
      af(np,4) = 0.5d0*betan(nps)/gamman(nps)         ! West
      af(np,5) = 1.d0 - 0.5d0*betan(nps)/gamman(nps)  ! P

      bf(np) = 0.d0

   end subroutine set_null_normal_grad_north

   !============================================================================

   !> Calculates the field over the north boundary based on the null normal gradient boundary condition
   subroutine set_null_normal_grad_north_5d(nx, ny, gamman, betan, af, bf)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: af(nx*ny,5)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np, nps

      j = ny

      iaux = nx*(j - 1)

      ! Calculating coefficients for the fictitious volume of the corner NW

      i = 2

      np = iaux + i
      nps = np - nx

      af(np,:) = 0.d0

      af(np,1) = -1.d0                          ! South
      af(np,3) = 1.d0 + betan(nps)/gamman(nps)  ! P
      af(np,4) = -betan(nps)/gamman(nps)        ! East

      bf(np) = 0.d0

      ! Calculating coefficients for the fictitious volumes of the north boundary, except NW and NE

      do i = 3, nx-2

         np = iaux + i
         nps = np - nx

         af(np,:) = 0.d0

         af(np,1) = -1.d0                          ! South
         af(np,2) = 0.5d0*betan(nps)/gamman(nps)   ! West
         af(np,3) = 1.d0                           ! P
         af(np,4) = -0.5d0*betan(nps)/gamman(nps)  ! East

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the fictitious volume of the corner NE

      i = nx - 1

      np = iaux + i
      nps = np - nx

      af(np,:) = 0.d0

      af(np,1) = -1.d0                          ! South
      af(np,2) = betan(nps)/gamman(nps)         ! West
      af(np,3) = 1.d0 - betan(nps)/gamman(nps)  ! P

      bf(np) = 0.d0

   end subroutine set_null_normal_grad_north_5d

   !============================================================================

   !> Calculates the field over the south boundary based on the null normal gradient boundary condition
   subroutine set_null_normal_grad_south(nx, ny, gamman, betan, af, bf)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np

      j = 1

      iaux = nx*(j - 1)

      ! Calculating coefficients for the fictitious volume of the corner SW

      i = 2

      np = iaux + i

      af(np,:) = 0.d0

      af(np,5) = -1.d0 + 0.5d0*betan(np)/gamman(np)   ! P
      af(np,6) = -0.5d0*betan(np)/gamman(np)          ! East
      af(np,8) = 1.d0 + 0.5d0*betan(np)/gamman(np)    ! North
      af(np,9) = -0.5d0*betan(np)/gamman(np)          ! Northeast

      bf(np) = 0.d0

      ! Calculating coefficients for the fictitious volumes of the south boundary, except SW and SE

      do i = 3, nx-2

         np = iaux + i

         af(np,:) = 0.d0

         af(np,4) = 0.25d0*betan(np)/gamman(np)    ! West
         af(np,5) = -1.d0                          ! P
         af(np,6) = -0.25d0*betan(np)/gamman(np)   ! East
         af(np,7) = 0.25d0*betan(np)/gamman(np)    ! North-west
         af(np,8) = 1.d0                           ! North
         af(np,9) = -0.25d0*betan(np)/gamman(np)   ! Northeast

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the fictitious volume of the corner SE

      i = nx - 1

      np = iaux + i

      af(np,:) = 0.d0

      af(np,4) = 0.5d0*betan(np)/gamman(np)           ! West
      af(np,5) = -1.d0 - 0.5d0*betan(np)/gamman(np)   ! P
      af(np,7) = 0.5d0*betan(np)/gamman(np)           ! North-west
      af(np,8) = 1.d0 - 0.5d0*betan(np)/gamman(np)    ! North

      bf(np) = 0.d0

   end subroutine set_null_normal_grad_south

   !============================================================================

   !> Calculates the field over the south boundary based on the null normal gradient boundary condition
   subroutine set_null_normal_grad_south_5d(nx, ny, gamman, betan, af, bf)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: gamman(nx*ny)   !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(in) :: betan(nx*ny)    !< (metric tensor component) beta  at the center of north face of volume P (m2)
      real(8), intent(inout) :: af(nx*ny,5)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, iaux, np

      j = 1

      iaux = nx*(j - 1)

      ! Calculating coefficients for the fictitious volume of the corner SW

      i = 2

      np = iaux + i

      af(np,:) = 0.d0

      af(np,3) = -1.d0 + betan(np)/gamman(np)   ! P
      af(np,4) = -betan(np)/gamman(np)          ! East
      af(np,5) = 1.d0                           ! North

      bf(np) = 0.d0

      ! Calculating coefficients for the fictitious volumes of the south boundary, except SW and SE

      do i = 3, nx-2

         np = iaux + i

         af(np,:) = 0.d0

         af(np,2) = 0.5d0*betan(np)/gamman(np)  ! West
         af(np,3) = -1.d0                       ! P
         af(np,4) = -0.5d0*betan(np)/gamman(np) ! East
         af(np,5) = 1.d0                        ! North

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the fictitious volume of the corner SE

      i = nx - 1

      np = iaux + i

      af(np,:) = 0.d0

      af(np,2) = betan(np)/gamman(np)           ! West
      af(np,3) = -1.d0 - betan(np)/gamman(np)   ! P
      af(np,5) = 1.d0                           ! North

      bf(np) = 0.d0

   end subroutine set_null_normal_grad_south_5d

   !============================================================================

   !> Calculates the field over the west boundary based on the null normal gradient boundary condition
   subroutine set_null_normal_grad_west(nx, ny, alphae, betae, af, bf)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: alphae(nx*ny)   !< (metric tensor component) gamma at the center of east face of volume P (m2)
      real(8), intent(in) :: betae(nx*ny)    !< (metric tensor component) beta  at the center of east face of volume P (m2)
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, np

      i = 1

      ! Calculating coefficients for the fictitious volume of the corner SW

      j = 2

      np = nx*(j - 1) + i

      af(np,:) = 0.d0

      af(np,5) = -1.d0 + 0.5d0*betae(np)/alphae(np)     ! P
      af(np,6) = 1.d0 + 0.5d0*betae(np)/alphae(np)      ! East
      af(np,8) = -0.5d0*betae(np)/alphae(np)                  ! North
      af(np,9) = -0.5d0*betae(np)/alphae(np)                  ! Northeast

      bf(np) = 0.d0

      ! Calculating coefficients for the fictitious volumes of the west boundary, except SW and NW

      do j = 3, ny-2

         np = nx*(j - 1) + i

         af(np,:) = 0.d0

         af(np,2) = 0.25d0*betae(np)/alphae(np)    ! South
         af(np,3) = 0.25d0*betae(np)/alphae(np)    ! Southeast
         af(np,5) = -1.d0                          ! P
         af(np,6) = 1.d0                           ! East
         af(np,8) = -0.25d0*betae(np)/alphae(np)   ! North
         af(np,9) = -0.25d0*betae(np)/alphae(np)   ! Northeast

         bf(np) = 0.d0

      end do

      ! Calculating coefficients for the fictitious volume of the corner NW

      j = ny - 1

      np = nx*(j - 1) + i

      af(np,:) = 0.d0

      af(np,2) = 0.5d0*betae(np)/alphae(np)           ! South
      af(np,3) = 0.5d0*betae(np)/alphae(np)           ! Southeast
      af(np,5) = -1.d0 - 0.5d0*betae(np)/alphae(np)   ! P
      af(np,6) = 1.d0 - 0.5d0*betae(np)/alphae(np)    ! East

      bf(np) = 0.d0

   end subroutine set_null_normal_grad_west

   !============================================================================

   !> Calculates the coeficients for the corners
   subroutine set_coef_extrapolation_to_corners(nx, ny, f, af, bf)
      implicit none
      integer, intent(in) :: nx              !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny              !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: f(nx*ny)        !< (metric tensor component) gamma at the center of north face of volume P (m2)
      real(8), intent(inout) :: af(nx*ny,9)  !< Coefficients of the linear system for phi
      real(8), intent(inout) :: bf(nx*ny)    !< Source of the linear system for phi

      ! Inner variables
      integer :: i, j, np, npn, nps, npe, npw, npne, npnw, npse, npsw

      ! Corners (extrapolation)

      ! SW
      i = 1
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npe = np + 1
      npne = npn + 1
      af(np,:) = 0.d0
      af(np,5) = 1.d0
      bf(np) = (f(npn) + f(npe) + f(npne))*0.3333333333333333d0

      !SE
      i = nx
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npw = np - 1
      npnw = npn - 1
      af(np,:) = 0.d0
      af(np,5) = 1.d0
      bf(np) = (f(npn) + f(npw) + f(npnw))*0.3333333333333333d0

      ! NW
      i = 1
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npe = np + 1
      npse = nps + 1
      af(np,:) = 0.d0
      af(np,5) = 1.d0
      bf(np) = (f(nps) + f(npe) + f(npse))*0.3333333333333333d0

      ! NE
      i = nx
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npw = np - 1
      npsw = nps - 1
      af(np,:) = 0.d0
      af(np,5) = 1.d0
      bf(np) = (f(nps) + f(npw) + f(npsw))*0.3333333333333333d0

   end subroutine set_coef_extrapolation_to_corners

   !============================================================================

   !> \brief Calculates the contravariant velocities Uc and Vc over the real boundaries
   subroutine get_Ucbe_Vcbe(nx, ny, xe, ye, xke, yke, ube, vbe, Ucbe, Vcbe) ! Last two are output
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in the csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in the eta direction (real+fictitious)
      real(8), intent(in) :: xe(nx*ny)   !< x_eta at the center of east  face of volume P (m)
      real(8), intent(in) :: ye(nx*ny)   !< y_eta at the center of east  face of volume P (m)
      real(8), intent(in) :: xke(nx*ny)  !< x_csi at the center of east  face of volume P (m)
      real(8), intent(in) :: yke(nx*ny)  !< y_csi at the center of east  face of volume P (m)
      real(8), intent(in) :: ube(ny)  !< u over the faces of the  east boundary (m/s)
      real(8), intent(in) :: vbe(ny)  !< v over the faces of the  east boundary (m/s)
      real(8), intent(out) :: Ucbe(ny) !< Uc over the faces of the east  boundary (m2/s)
      real(8), intent(out) :: Vcbe(ny) !< Vc over the faces of the east  boundary (m2/s)

      ! Inner variables
      integer :: i, j, np

      ! East boundary

      i = nx - 1

      do j = 2, ny-1

         np = nx*(j - 1) + i

         Ucbe(j) = ube(j)*ye(np) - vbe(j)*xe(np)

         Vcbe(j) = vbe(j)*xke(np) - ube(j)*yke(np)

      end do

   end subroutine get_Ucbe_Vcbe

   !============================================================================

   subroutine get_uin_vin_pin_Tin_Mw(nx, ny, gamma, Rg, po, T0, u, uin, vin, &
         pin, Tin, Mw) ! Five last variables are output
      implicit none
      integer, intent(in) :: nx        !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny        !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: gamma     !< gamma = Cpo/Cvo in the chamber
      real(8), intent(in) :: Rg        !< Perfect gas constant
      real(8), intent(in) :: po        !< Stagnation pressure in the chamber
      real(8), intent(in) :: T0        !< Stagnation temperature in the chamber
      real(8), intent(in) :: u(nx*ny)  !< Cartesian velocity of the last iteraction
      real(8), intent(out) :: uin(ny)  !< Velocity u in the entrance
      real(8), intent(out) :: vin(ny)  !< Velocity v in the entrance
      real(8), intent(out) :: pin(ny)  !< Pressure in the entrance
      real(8), intent(out) :: Tin(ny)  !< Temperature in the entrance
      real(8), intent(out) :: Mw(ny)   !< Mach number in the entrance

      integer :: i, j, np, npe

      i = 1

      do j = 2, ny-1

         np = nx*(j - 1) + i

         npe = np + 1

         uin(j) = (u(np) + u(npe))*0.5d0

         vin(j) = 0.d0 ! From boundary conditions

         Tin(j) = T0 - (gamma-1.d0)/(2.d0*gamma*Rg)*(uin(j)**2 + vin(j)**2)

         Mw(j) = dsqrt((uin(j)**2 + vin(j)**2)/(gamma*Rg*Tin(j)))

         pin(j) = po*(1.d0 + 0.5d0*(gamma-1.d0)*Mw(j)**2)**(-gamma/(gamma-1.d0))

      end do

   end subroutine get_uin_vin_pin_Tin_Mw

   !============================================================================

   subroutine get_plin_and_p_fictitious(nx, ny, pina, pin, plin, p) ! Output: last two entries
      implicit none
      integer, intent(in) :: nx        !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny        !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: pina(ny)  !< Pressure in the entrance at a time step before
      real(8), intent(in) :: pin(ny)   !< Pressure in the entrance
      real(8), intent(out) :: plin(ny) !< Pressure correction in the entrance
      real(8), intent(out) :: p(nx*ny) !< Pressure at center of volumes

      integer :: i, j, np, npe

      ! West boundary

      i = 1

      do j = 2, ny-1

         np = nx*(j - 1) + i

         npe = np + 1

         p(np) = 2.d0*pin(j) - p(npe)

         plin(j) = pin(j) - pina(j)

      end do

   end subroutine get_plin_and_p_fictitious

   !============================================================================

   subroutine get_u_v_extrapolation_to_fictitious(nx, ny, modvis, u, v) ! InOutput: last two entries
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in csi direction (real+fictitious)
      integer, intent(in) :: ny     !< Number of volumes in eta direction (real+fictitious)
      integer, intent(in) :: modvis !< Viscosity model (0=Euler; 1=Navier-Stokes)
      real(8), intent(inout) :: u(nx*ny)  !< Cartesian velocity
      real(8), intent(inout) :: v(nx*ny)  !< Cartesian velocity

      integer :: i, j, np, npe, npw, nps, npn, iaux

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         u(np) = u(npe)
         v(np) = -v(npe)
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         u(np) = u(npw)
         v(np) = v(npw)
      end do

      ! South boundary
      j = 1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         npn = np + nx
         u(np) = u(npn)
         v(np) = -v(npn)
      end do

      ! North boundary
      j = ny
      iaux = nx*(j - 1)
      if (modvis == 1) then ! Viscous

         do i = 2, nx-1
            np = iaux + i
            nps = np - nx
            u(np) = -u(nps)
            v(np) = -v(nps)
         end do

      else ! Inviscid

         do i = 2, nx-1
            np = iaux + i
            nps = np - nx
            u(np) = u(nps)
            v(np) = v(nps)
         end do

      end if

      ! Corners extrapolation

      call get_extrapolation_to_corners(nx, ny, u)

      call get_extrapolation_to_corners(nx, ny, v)

   end subroutine get_u_v_extrapolation_to_fictitious

   !============================================================================

   !> \brief Extrapolation of a field to the fictitious volumes of the corners
   subroutine get_extrapolation_to_corners(nx, ny, f)
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in the csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in the eta direction (real + fictitious)
      real(8), intent(inout) :: f(nx*ny) !< Field to be extrapolated to the corners

      integer :: i, j, np, npe, npw, npn, nps, npne, npnw, npse, npsw

      ! Corners (extrapolation)

      ! SW
      i = 1
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npe = np + 1
      npne = npn + 1
      f(np) = (f(npn) + f(npe) + f(npne))*0.3333333333333333d0

      !SE
      i = nx
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npw = np - 1
      npnw = npn - 1
      f(np) = (f(npn) + f(npw) + f(npnw))*0.3333333333333333d0

      ! NW
      i = 1
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npe = np + 1
      npse = nps + 1
      f(np) = (f(nps) + f(npe) + f(npse))*0.3333333333333333d0

      ! NE
      i = nx
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npw = np - 1
      npsw = nps - 1
      f(np) = (f(nps) + f(npw) + f(npsw))*0.3333333333333333d0

   end subroutine get_extrapolation_to_corners

   !============================================================================

   subroutine get_boundary_simplec_coefficients(nx, ny, de, dn) ! InOutput: de, dn
      implicit none
      integer, intent(in) :: nx  !< number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< number of volumes in eta direction (real + fictitious)
      real(8), intent(inout) :: de(nx*ny) !< SIMPLEC coefficient for Uce
      real(8), intent(inout) :: dn(nx*ny) !< SIMPLEC coefficient for Vcn

      integer :: i, j, np, npe, npw, iaux

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         de(np) = de(npe)
      end do

      ! East boundary
      i = nx-1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         de(np) = de(npw)
      end do

      ! South boundary
      j = 1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         dn(np) = 0.d0
      end do

      ! North boundary
      j = ny-1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         dn(np) = 0.d0
      end do

   end subroutine get_boundary_simplec_coefficients

   !============================================================================

   subroutine get_Uce_Vcn_at_boundary_faces(nx, ny, ye, u, Uce, Vcn) ! Two last are output
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: ye(nx*ny)    !< y_eta at face east of volume P
      real(8), intent(in) :: u(nx*ny)     !< Cartesian velocity of the last iteraction
      real(8), intent(out) :: Uce(nx*ny)  !< Contravariant velocity U at east face
      real(8), intent(out) :: Vcn(nx*ny)  !< Contravariant velocity V at north face

      integer :: i, j, np, npe, iaux

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         Uce(np) = 0.5d0*(u(npe) + u(np))*ye(np)
      end do

      ! East boundary
      i = nx-1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         Uce(np) = 0.5d0*(u(npe) + u(np))*ye(np)
      end do

      ! South boundary
      j = 1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         Vcn(np) = 0.d0
      end do

      ! North boundary
      j = ny-1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np   = iaux + i
         Vcn(np) = 0.d0
      end do

   end subroutine get_Uce_Vcn_at_boundary_faces

   !============================================================================

   subroutine get_Uce_Vcn_at_boundary_faces_with_pl(nx, ny, pl, de, Uce) ! InOutput: Last entry
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(in) :: pl(nx*ny)       !< Pressure correction
      real(8), intent(in) :: de(nx*ny)       !< SIMPLEC coefficient for Uce
      real(8), intent(inout) :: Uce(nx*ny)   !< Contravariant velocity U at east face

      integer :: i, j, np, npe

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         Uce(np) = Uce(np) + de(np)*(pl(np) - pl(npe))
      end do

      ! West boundary
      i = nx-1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         Uce(np) = Uce(np) + de(np)*(pl(np) - pl(npe))
      end do

      ! North and South Vcn do not neet to be changed

   end subroutine get_Uce_Vcn_at_boundary_faces_with_pl

   !============================================================================

   subroutine get_isentropic_mass_flow(po, T0, gamma, Rg, Sg, fm1D)
      implicit none
      real(8), intent(in) :: po     !< Stagnation pressure
      real(8), intent(in) :: T0     !< Stagnation temperature
      real(8), intent(in) :: gamma  !< Specific heat ratio
      real(8), intent(in) :: Rg     !< Perfect gas constant
      real(8), intent(in) :: Sg     !< Throttle area
      real(8), intent(out) :: fm1D  !< Mass flow rate

      fm1D = po*Sg*sqrt(gamma/(Rg*T0)*(2.d0/(gamma + 1.d0))**((gamma &
         + 1.d0)/(gamma - 1.d0)))

   end subroutine get_isentropic_mass_flow

   !============================================================================

   subroutine get_mach_area(kf, ar, gamma, M)
      implicit none
      integer, intent(in) :: kf     !< Kind of flow (0=subsonic; 1=supersonic)
      real(8), intent(in) :: ar     !< Area ratio (local area/throttle area)
      real(8), intent(in) :: gamma  !< Specific heat ratio
      real(8), intent(out) :: M     !< Mach number at the local area

      integer :: it
      real(8) :: M0, M1, f, x

      f(x) = (2.d0*(1.d0 + (gamma - 1.d0)*x**2/2.d0)/(gamma + 1.d0))**((gamma &
         + 1.d0)/(gamma - 1.d0)) - (x*ar)**2

      if (abs(ar-1.d0) < 1.d-15) then ! Sonic flow in the throttle

         M = 1.d0

      else

         if (kf == 0) then ! Subsonic
            M0 = 1.d-4
            M1 = 1.d0
         else              ! Supersonic
            M0 = 1.d0
            M1 = 100.d0
         end if

         do it = 1, 1000

            M = 0.5d0*(M0 + M1)

            if (f(M)*f(M0) > 0) then
               M0 = M
            else
               M1 = M
            end if

            if (abs(M1-M0) < 1.d-15) exit

         end do

      end if

   end subroutine get_mach_area

   !============================================================================

   subroutine get_initial_guess(nx, ny, ig, modvis, beta, po, T0, gamma, Rg, &
         Sg, ye, yk, radius, rn, & ! Input
      M1D, p1D, T1D, u1D, p, T, u, v, ue, un, Uce, Vcn, uin, vin, pin, Tin, &
         Mw, fm1D, Fd1D, Fpv1D, de, dn, ro, roe, ron, a, ap, b) ! Output
      implicit none
      integer, intent(in) :: nx           !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny           !< Number of volumes in eta direction (real + fictitious)
      integer, intent(in) :: ig           !< Volume number that has to face east on the throat of the nozzle
      integer, intent(in) :: modvis       !< Viscosity model (0=Euler; 1=Navier-Stokes)
      real(8), intent(in) :: beta         !< Constant of the UDS/CDS mixing scheme
      real(8), intent(in) :: po           !< Stagnation pressure in the chamber
      real(8), intent(in) :: T0           !< Stagnation temperature in the chamber
      real(8), intent(in) :: gamma        !< gamma = Cpo/Cvo in the chamber
      real(8), intent(in) :: Rg           !< Perfect gas constant
      real(8), intent(in) :: Sg           !< Throttle area
      real(8), intent(in) :: ye(nx*ny)    !< y_eta at face east of volume P
      real(8), intent(in) :: yk(nx*ny)    !< y_csi at face north of volume P
      real(8), intent(in) :: radius(nx*ny)!< Radius of northest corner of volume P
      real(8), intent(in) :: rn(nx*ny)    !< Radius of the center of north face of volume P
      real(8), intent(out) :: fm1D
      real(8), intent(out) :: Fd1D
      real(8), intent(out) :: Fpv1D
      real(8), intent(out) :: M1D(nx)     !< Mach number of the isentropic flow
      real(8), intent(out) :: p1D(nx)     !< Pressure of the isentropic flow
      real(8), intent(out) :: T1D(nx)     !< Temperature of the isentropic flow
      real(8), intent(out) :: u1D(nx)     !< Velocity of the isentropic flow
      real(8), intent(out) :: p(nx*ny)    !< Pressure at center o volume P
      real(8), intent(out) :: T(nx*ny)    !< Temperature of the last iteraction
      real(8), intent(out) :: u(nx*ny)    !< Cartesian velocity of the last iteraction
      real(8), intent(out) :: v(nx*ny)    !< Cartesian velocity of the last iteraction
      real(8), intent(out) :: ue(nx*ny)   !< Cartesian velocity u at center of east face
      real(8), intent(out) :: un(nx*ny)   !< Cartesian velocity u at center of north face
      real(8), intent(out) :: Uce(nx*ny)  !< Contravariant velocity U at east face
      real(8), intent(out) :: Vcn(nx*ny)  !< Contravariant velocity V at north face
      real(8), intent(out) :: uin(ny)     !< Velocity u in the entrance
      real(8), intent(out) :: vin(ny)     !< Velocity v in the entrance
      real(8), intent(out) :: pin(ny)     !< Pressure in the entrance
      real(8), intent(out) :: Tin(ny)     !< Temperature in the entrance
      real(8), intent(out) :: Mw(ny)      !< Mach number in the entrance
      real(8), intent(out) :: de(nx*ny)   !< Simplec coef. for the contravariant velocity U (east face)
      real(8), intent(out) :: dn(nx*ny)   !< Simplec coef. for the contravariant velocity V (north face)
      real(8), intent(out) :: ro(nx*ny)   !< Specific mass (absolute density) at center of volumes
      real(8), intent(out) :: roe(nx*ny)  !< Absolute density at east face
      real(8), intent(out) :: ron(nx*ny)  !< Absolute density at north face
      real(8), intent(out) :: a(nx*ny,9)  !< Coefficients of the linear system for u, v, T
      real(8), intent(out) :: ap(nx*ny,5) !< Coefficients of the linear system
      real(8), intent(out) :: b(nx*ny)    !< Source vector of the linear system

      integer :: i, j, np, npe, npw, kf
      real(8), parameter :: pi = dacos(-1.d0)
      real(8) :: Si, Me1D, pe1D, Te1D, ue1D, aux

      call get_isentropic_mass_flow(po, T0, gamma, Rg, Sg, fm1D)

      ! Calculating isentropic solution
      do i = 2, nx-1

         j = ny-1

         np = nx*(j - 1) + i

         Si = pi*rn(np)**2

         kf = 0

         if (i > ig) kf = 1

         call get_mach_area(kf, Si/Sg, gamma, M1D(i))

         aux = 1.d0 + 0.5d0*(gamma-1.d0)*M1D(i)**2

         p1D(i) = po*aux**(-gamma/(gamma-1.d0))

         T1D(i) = T0/aux

         u1D(i) = M1D(i)*sqrt(gamma*Rg*T1D(i))

         do j = 1, ny

            np = nx*(j - 1) + i

            p(np) = p1D(i)

            T(np) = T1D(i)

            u(np) = u1D(i)

            if (j > 1 .and. j < ny-1) then
               un(np) = u1D(i)
               Vcn(np) = -un(np)*yk(np) ! vn = 0
            end if

         end do

      end do

      ! Calculates u and v at fictitious nodes using boundary conditions
      call get_u_v_extrapolation_to_fictitious(nx, ny, modvis, u, v) ! InOutput: last two entries

      ! Calculates gas proprieties in the entrance
      call get_uin_vin_pin_Tin_Mw(nx, ny, gamma, Rg, po, T0, u & ! Input
      , uin, vin, pin, Tin, Mw)  ! Output

      ! Calculates p and T in the fictitious nodes of the entrance according to boundary conditions
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         p(np) = 2.d0*pin(j) - p(npe)
         T(np) = 2.d0*Tin(j) - T(npe)
      end do

      ! Calculates p and T in the fictitious nodes of the exit according to boundary conditions
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         p(np) = p(npw)
         T(np) = T(npw)
      end do

      ! Initial guess of velocities in east face of each CV
      ! and analytical solution in the entrance and exit boundaries
      do i = 1, nx-1

         j = ny-1

         np = nx*(j - 1) + i

         Si = pi*radius(np)**2

         kf = 0

         if (i > ig) kf = 1

         call get_mach_area(kf, Si/Sg, gamma, Me1D)

         aux = 1.d0 + 0.5d0*(gamma-1.d0)*Me1D**2

         pe1D = po*aux**(-gamma/(gamma-1.d0))

         Te1D = T0/aux

         ue1D = Me1D*sqrt(gamma*Rg*Te1D)

         do j = 2, ny-1
            np = nx*(j - 1) + i
            ue(np) = ue1D
            Uce(np) = ue1D*ye(np) ! v = 0
         end do

         ! Entrance boundary
         if (i == 1) then
            u1D(i) = ue1D
            p1D(i) = pe1D
            T1D(i) = Te1D
            M1D(i) = Me1D
         end if

         ! Exit boundary
         if (i == nx-1) then
            u1D(nx) = ue1D
            p1D(nx) = pe1D
            T1D(nx) = Te1D
            M1D(nx) = Me1D

            Fd1D = fm1D*ue1D
            Fpv1D = pe1D*Si
         end if

      end do

      ! Contravariant velocities at boundary faces
      call get_Uce_Vcn_at_boundary_faces(nx, ny, ye, u, Uce, Vcn) ! Output: Uce, Vcn

      ! Simplec coefficients at boundary faces
      call get_boundary_simplec_coefficients(nx, ny, de, dn) ! InOutput: de, dn

      ! Specific mass at nodes
      call get_density_at_nodes(nx, ny, Rg, p, T, ro) ! ro is output

      ! Specific mass at faces
      call get_density_at_faces(nx, ny, beta, ro, Uce, Vcn, roe, ron) ! roe and ron are output

      a = 0.d0
      ap = 0.d0
      b = 0.d0

      ap(:,3) = 1.d0

   end subroutine get_initial_guess

   ! ===========================================================================

   subroutine get_boundary_nodes(folder_input, folder_output, sim_id, nmbr, &
         optm, lid, reload, nx, ctd, ig, sg, rcgd, rg) ! Output: last three entries
      implicit none
      character(3), intent(out) :: nmbr               !< process name in parallel processing
      character(200), intent(in) :: folder_input      !< input folder name
      character(200), intent(inout) :: folder_output  !< output folder name
      character(200), intent(inout) :: sim_id         !< simulation name
      integer, intent(in) :: reload
      integer, intent(in) :: optm   !< Coupling with optimization method (0=no; 1=yes)
      integer, intent(in) :: lid    !< Listing file id
      integer, intent(in) :: nx     !< Number of volumes in the axial direction
      integer, intent(in) :: ctd    !< Kind of grid in axial direction (0=uniform; 1=concentrated)
      integer, intent(out) :: ig    !< Volume number that has to face east on the throat of the nozzle
      real(8), intent(out) :: sg    !< Area of the throttle
      real(8), intent(out) :: rcgd  !< Radius of curvature in the throat in the divergent side
      real(8), intent(out) :: rg    !< Radius in the throat of the nozzle [m]

      character(10), dimension(:), allocatable :: xname  !< name of variables
      character(200) :: sname   !< name of the simulation
      integer :: ind            !< number of the individual
      integer, dimension(:), allocatable :: xopt         !< checker of variables to be read from DEPP

      integer :: i
      integer :: k
      integer :: prt    !< Print nozzle coefficients (0=no; 1=yes)
      integer :: cont   !< count the number of parameters read
      integer :: ns     !< Nozzle shape (1=arbitrary; 2=conical; 3=parabolic; 4=bell-arc/spline; 5=bell-spline)
      integer :: npar   !< Number of the parameters
      integer :: ncp    !< Number of control points
      real(8) :: dist   !< distance between two points
      real(8) :: raio   !< local radius
      real(8) :: area   !< area of the surface wall (prescribed only if greater than zero)
      real(8) :: rl     !< Local radius [m]
      real(8) :: dx     !< Length of volume [m]
      real(8) :: lcc    !< Length of the chamber + convergente section [m]
      real(8) :: x0     !< Coordinated at the entrance of the chamber [m]
      real(8) :: x1     !< Coordinated in input of the region convergent [m]
      real(8) :: ln     !< Length of the nozzle (convergent + divergent) [m]
      real(8) :: l      !< Total length (chamber + nozzle) [m]
      real(8) :: ri     !< Radius in the entrance of the nozzle [m]
      real(8) :: x(nx)  !< x coordinates of the north contour of the nozzle [m]
      real(8) :: y(nx)  !< y coordinates of the north contour of the nozzle [m]
      real(8) :: incl   !< Inclination of wall of the convergent region [degrees]
      real(8) :: rcin   !< Curvature radius in input [m]
      real(8) :: rcgc   !< Radius of curvature in the throat in the convergent side
      real(8) :: lch    !< Length of the chamber [m]
      real(8) :: lc     !< Length of the convergente section [m]
      real(8) :: cg     !< Coefficient of concentration of the mesh
      real(8) :: xc4    !< x coordinate of arc center of the nozzle throat [m]
      real(8) :: yc4d   !< y coordinate of the center of circumference of the curvature of throat in the divergent side
      real(8) :: x4     !< x coordinate of the nozzle throat [m]
      real(8) :: y4     !< y coordinate of the nozzle throat [m]
      real(8) :: x5     !< x coordinate in the intersection between the arc and spline [m]
      real(8) :: y5     !< y coordinate in the intersection between the arc and spline [m]
      real(8) :: dy5    !< Inclination of wall in the intersection between the arc and spline
      real(8) :: xaux
      real(8) :: yaux
      real(8), dimension(:), allocatable :: par    !< Parameters that define the nozzle (length, inclination and control points)
      real(8), dimension(:), allocatable :: xj     !< x coordinates of the control points
      real(8), dimension(:,:), allocatable :: spl  !< Coefficients of the spline
      real(8), parameter :: pi = dacos(-1.d0)

      logical :: lexist
      integer :: dos

      open(49, file = trim(folder_input) // 'geometry.dat')

      read(49,*) ns
      read(49,*) npar
      if (ns == 2) npar = 2
      if (ns == 3) npar = 3
      if (ns == 6) npar = 0
      read(49,*) incl
      read(49,*) rcin
      read(49,*) rcgc
      read(49,*) rcgd
      read(49,*) rg
      read(49,*) area
      read(49,*) lch
      read(49,*) lc

      ncp = 0

      prt = 0

      x0 = -lch - lc
      x1 = -lc
      lcc = abs(x0)

      if (npar > 0) allocate(par(npar))
      if (npar > 0) allocate(xopt(npar))
      if (npar > 0) allocate(xname(npar))

      cont = 0

      do i = 1, npar

         if (optm /= 0) then

            read(49,*) par(i), xopt(i), xname(i)

         else

            read(49,*) par(i)

         end if

         cont = cont + 1

      end do

      ind = 0

      call convert_int_to_char3(ind, nmbr)

      if (optm /= 0) then

         call depp_get_parameters(npar, xopt, xname, par, ind, sname)

         call convert_int_to_char3(ind, nmbr)

         sim_id = trim(adjustl(sim_id)) // nmbr

      end if

      if (ns == 1) then ! nozzle arbitrarily defined by arbitrary_contour

         if (area > 0) call get_length_of_the_nozzle(lid, prt, optm, ns, npar, &
            ncp, rg, rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)

         ln = par(1) + dabs(x1)
         l = dabs(x0) + par(1)

         x(1) = x0
         y(1) = arbitrary_contour(x0, x1, x(1), incl, rcin, rcgc, rcgd, rg, par)

         x(nx-1) = par(1)
         y(nx-1) = arbitrary_contour(x0, x1, x(nx-1), incl, rcin, rcgc, rcgd, &
            rg, par)

         dx = (x(nx-1) - x(1))/(nx - 2)

         do i = 2, nx-2
            x(i) = x(1) + (i - 1)*dx
            y(i) = arbitrary_contour(x0, x1, x(i), incl, rcin, rcgc, rcgd, rg, &
               par)
         end do

         x4 = 0.d0
         y4 = arbitrary_contour(x0, x1, x4, incl, rcin, rcgc, rcgd, rg, par)

         xaux = x4 + 1.d-8
         yaux = arbitrary_contour(x0, x1, xaux, incl, rcin, rcgc, rcgd, rg, par)

         rg = y4
         rcgd = ((xaux - x4)**2 + (yaux - y4)**2)/(2*(yaux - y4))

      else if (ns == 2) then ! conic nozzle

         par(2) = par(2)*pi/180.d0

         if (area > 0) call get_length_of_the_nozzle(lid, prt, optm, ns, npar, &
            ncp, rg, rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)

         ln = par(1) + dabs(x1)
         l = dabs(x0) + par(1)

         x(1) = x0
         y(1) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, x(1))

         x(nx-1) = par(1)
         y(nx-1) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, x(nx-1))

         dx = (x(nx-1) - x(1))/(nx - 2)

         do i = 2, nx-2
            x(i) = x(1) + (i - 1)*dx
            y(i) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, x(i))
         end do

      else if (ns == 3) then ! parabolic nozzle

         if (area > 0) call get_length_of_the_nozzle(lid, prt, optm, ns, npar, &
            ncp, rg, rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)

         ln = par(1) + dabs(x1)
         l = dabs(x0) + par(1)

         x(1) = x0
         y(1) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, x1, &
            incl, x(1))

         x(nx-1) = par(1)
         y(nx-1) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, x1, &
            incl, x(nx-1))

         dx = (x(nx-1) - x(1))/(nx - 2)

         do i = 2, nx-2
            x(i) = x(1) + (i - 1)*dx
            y(i) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, x1, &
               incl, x(i))
         end do

      else if (ns == 4) then ! bell nozzle with arc and spline

         ncp = npar - 2

         allocate(xj(ncp), spl(ncp-1,4))

         xj(1) = par(2)
         if (ncp > 2) xj(2) = 0.1d0*par(1)
         if (ncp > 3) xj(3) = 0.3d0*par(1)
         if (ncp > 4) xj(4) = 0.6d0*par(1)
         xj(ncp) = par(1)

         xc4 = 0.d0
         yc4d = rcgd + rg

         ! intersection between the arc and spline
         x5 = par(2)
         y5 = -sqrt(rcgd**2 - (x5 - xc4)**2) + yc4d
         dy5 = (x5 - xc4)/sqrt(rcgd**2 - (x5 - xc4)**2)

         call get_spline_coefficients(y5, dy5, par(npar), xj, par(3:npar-1), &
            ncp-1, spl)

         if (area > 0) call get_length_of_the_nozzle(lid, prt, optm, ns, npar, &
            ncp, rg, rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)

         ln = par(1) + dabs(x1)
         l = dabs(x0) + par(1)

         x(1) = x0
         y(1) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
            spl, ncp+2, xj, x(1))

         x(nx-1) = par(1)
         y(nx-1) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
            spl, ncp+2, xj, x(nx-1))

         dx = (x(nx-1) - x(1))/(nx - 2)

         do i = 2, nx-2
            x(i) = x(1) + (i - 1)*dx
            y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
               spl, ncp+2, xj, x(i))
         end do

      else if (ns == 5) then ! bell nozzle with spline function

         ncp = npar - 1

         allocate(xj(ncp), spl(ncp-1,4))

         xj(1) = 0.d0
         xj(2) = 0.04d0*par(1)
         if (ncp > 3) xj(3) = 0.1d0*par(1)
         if (ncp > 4) xj(4) = 0.3d0*par(1)
         if (ncp > 5) xj(5) = 0.6d0*par(1)
         xj(ncp) = par(1)

         x5 = 0.d0

         call get_spline_coefficients(rg, dy5, par(npar), xj, par(2:ncp), &
            ncp-1, spl)

         if (area > 0) call get_length_of_the_nozzle(lid, prt, optm, ns, npar, &
            ncp, rg, rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)

         ln = par(1) + dabs(x1)
         l = dabs(x0) + par(1)

         x(1) = x0
         y(1) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
            spl, ncp+2, xj, x(1))

         x(nx-1) = par(1)
         y(nx-1) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
            spl, ncp+2, xj, x(nx-1))

         dx = (x(nx-1) - x(1))/(nx - 2)

         do i = 2, nx-2
            x(i) = x(1) + (i - 1)*dx
            y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
               spl, ncp+2, xj, x(i))
         end do

      else if (ns == 6) then ! nozzle defined by the geometry_points.dat file

         open(87, file = trim(folder_input) // 'geometry_points.dat')

         read(87,*) ncp

         npar = ncp + 1

         allocate(xj(ncp), spl(ncp-1,4), par(npar))

         read(87,*) xj(1), yaux

         if (dabs(rg - yaux) > 1.d-8) then

            write(*,*) "= ERROR: Radius informed in the files 'geometry.dat' ="
            write(*,*) "====== and 'geometry_points.dat' are differents. ====="

            stop

         end if

         do i = 2, ncp

            read(87,*) xj(i), par(i)

         end do

         par(1) = xj(ncp)

         par(npar) = 200.d0

         x5 = 0.d0

         call get_spline_coefficients(rg, dy5, par(npar), xj, par(2:ncp), &
            ncp-1, spl)

         if (area > 0) call get_length_of_the_nozzle(lid, prt, optm, ns, npar, &
            ncp, rg, rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)

         ln = par(1) + dabs(x1)
         l = dabs(x0) + par(1)

         x(1) = x0
         y(1) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
            spl, ncp+2, xj, x(1))

         x(nx-1) = par(1)
         y(nx-1) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
            spl, ncp+2, xj, x(nx-1))

         dx = (x(nx-1) - x(1))/(nx - 2)

         do i = 2, nx-2
            x(i) = x(1) + (i - 1)*dx
            y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
               spl, ncp+2, xj, x(i))
         end do

         !         call convert_int_to_char3(ncp-1, nmbr)
         !         open(98, file = 'geometry_points_t2000_' // nmbr // '.dat')
         !         write(98,*) 1.d0/(ncp-1)
         !         dx = 0.025d0
         !         do i = 1, 40
         !            yaux = i*dx
         !            yaux = bell_contour(lid, prt, rg, rcg, rcin, x1, x5, incl, spl, &
            !               ncp+2, xj, yaux)
         !            write(98,*) i, yaux, i*dx
         !         end do

      else

         write(*,*) " ======================================================== "
         write(*,*) " ====== ERROR: (NS) INVALID OPTION FOR NOZZLE SHAPE ===== "
         write(*,*) " ====== 1=ARBITRARY; 2=CONIC; 3=PARABOLIC; 4:5=BELL ===== "
         write(*,*) " ======================================================== "
         stop

      end if

      ! concentration of the volumes near throat
      if (ctd/= 0) then

         do k = 1, 100

            cg = l/(0.5d0*y(1) + sum(y(2:nx-2)) + 0.5d0*y(nx-1))

            do i = 2, nx-2

               dx = cg*(y(i-1) + y(i))/2.d0
               x(i) = x(i-1) + dx

               select case (ns)
                  case (1)
                     y(i) = arbitrary_contour(x0, x1, x(i), incl, rcin, rcgc, &
                        rcgd, rg, par)
                  case (2)
                     y(i) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, &
                        incl, x(i))
                  case (3)
                     y(i) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, &
                        rcin, x1, incl, x(i))
                  case(4,5,6)
                     y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, &
                        x5, incl, spl, ncp+2, xj, x(i))
                  case default
                     write(*,*) "ERROR(1): invalid value for NS."
               end select

            end do

         end do

      end if

      ! identification of the volume that has the east side nearest the throat
      rl = y(1)

      do i = 1, nx-1

         if (y(i) < rl) then
            rl = y(i)
            ig = i
         end if

      end do

      ! adjustment to the eastern side of the volume IG is on top of the throat
      if (ctd == 0) then

         x(ig) = 0.d0
         dx = lcc/(ig - 1.d0)

         do i = 2, ig

            x(i) = x(1) + (i - 1.d0)*dx

            select case (ns)
               case (1)
                  y(i) = arbitrary_contour(x0, x1, x(i), incl, rcin, rcgc, &
                     rcgd, rg, par)
               case (2)
                  y(i) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, &
                     x(i))
               case (3)
                  y(i) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, &
                     rcin, x1, incl, x(i))
               case (4,5,6)
                  y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, &
                     incl, spl, ncp+2, xj, x(i))
               case default
                  write(*,*) "ERROR(4): invalid value for NS."
            end select

         end do

         dx = (l - lcc)/(nx - ig - 1.d0)

         yaux = 0.d0

         do i = ig+1, nx-2

            x(i) = x(ig) + (i - ig)*dx

            select case (ns)
               case (1)
                  y(i) = arbitrary_contour(x0, x1, x(i), incl, rcin, rcgc, &
                     rcgd, rg, par)
               case (2)
                  y(i) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, &
                     x(i))
               case (3)
                  y(i) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, &
                     rcin, x1, incl, x(i))
               case (4,5,6)
                  y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, &
                     incl, spl, ncp+2, xj, x(i))
               case default
                  write(*,*) "ERROR(5): invalid value for NS."
            end select

            if (optm /= 0 .and. y(i) < yaux) then

               call depp_save_fitness(0.d0, 2, "ERROR: convergent wall")

               stop

            end if

            yaux = y(i)

         end do

      else

         x(ig) = 0.d0

         do k = 1, 100

            cg = lcc/(0.5d0*y(1) + sum(y(2:ig-1)) + 0.5d0*y(ig))

            do i = 2, ig

               dx = cg*(y(i-1) + y(i))/2.d0
               x(i) = x(i-1) + dx

               select case (ns)
                  case (1)
                     y(i) = arbitrary_contour(x0, x1, x(i), incl, rcin, rcgc, &
                        rcgd, rg, par)
                  case (2)
                     y(i) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, &
                        incl, x(i))
                  case (3)
                     y(i) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, &
                        rcin, x1, incl, x(i))
                  case (4,5,6)
                     y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, &
                        x5, incl, spl, ncp+2, xj, x(i))
                  case default
                     write(*,*) "ERROR(2): invalid value for NS."
               end select

            end do

         end do

         do k = 1, 100

            yaux = 0.d0

            cg = (l-lcc)/(0.5*y(ig) + sum(y(ig+1:nx-2)) + 0.5*y(nx-1))

            do i = ig+1, nx-2

               dx = cg*(y(i-1) + y(i))/2.d0
               x(i) = x(i-1) + dx

               select case (ns)
                  case (1)
                     y(i) = arbitrary_contour(x0, x1, x(i), incl, rcin, rcgc, &
                        rcgd, rg, par)
                  case (2)
                     y(i) = conical_contour(par, rg, rcgc, rcgd, rcin, x1, &
                        incl, x(i))
                  case (3)
                     y(i) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, &
                        rcin, x1, incl, x(i))
                  case (4,5,6)
                     y(i) = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, &
                        x5, incl, spl, ncp+2, xj, x(i))
                  case default
                     write(*,*) "ERROR(3): invalid value for NS."
               end select

               if (optm /= 0 .and. y(i) < yaux) then

                  call depp_save_fitness(0.d0, 2, "ERROR: convergent wall")

                  stop

               end if

               yaux = y(i)

            end do

         end do

      end if

      ! calculates the wall surface
      area = 0.d0
      do i = ig, nx-2
         dist = dsqrt((x(i+1) - x(i))**2 + (y(i+1) - y(i))**2)
         raio = 0.5d0*(y(i) + y(i+1))
         area = area + 2.d0*pi*raio*dist
      end do

      ri = y(1)
      sg = pi*rg**2

      ! Openning listing file

      if (optm /= 0) then

         folder_output = trim(folder_output) // trim(adjustl(sname)) // "/"

      end if

      folder_output = trim(folder_output) // trim(adjustl(sim_id)) // "/"

      inquire(file = trim(folder_output), exist = lexist)

      if (.not. lexist) then

         dos = system("mkdir " // trim(folder_output))

      else

         if (reload == 0) then

            dos = system("rm -r " // trim(folder_output))

            dos = system("mkdir " // trim(folder_output))

         end if

      end if

      inquire(file = trim(folder_output) // trim(adjustl(sim_id)) // &
         ".lst", exist = lexist)

      if (lexist .and. reload == 1) then
         open(lid, file = trim(folder_output) // &
            trim(adjustl(sim_id)) // ".lst", POSITION = 'APPEND')

         write(lid,*)
         write(lid,*) " ====================================================", &
            "========================== "
         write(lid,*)
         write(lid,*) "           RELOADED. CONTINUING THE CALCULATIONS..."
         write(lid,*)
         write(lid,*) " ====================================================", &
            "========================== "
         write(lid,*)

      else
         open(lid, file = trim(folder_output) // trim(adjustl(sim_id)) // &
            ".lst")
      end if

      write(lid,*)
      write(lid,*) "LISTING FILE OF MACH2D"
      write(lid,*)

      ! Write geometry parameters

      write(lid,*)
      write(lid,*) " ==========================  GEOMETRY PARAMETERS  ", &
         "==========================="
      write(lid,*)
      write(lid,"(i15, a, a)") ns, " = ns:     Nozzle shape (1=arbitrary; ", &
         "2=conical; 3=parabolic; 4=bell-arc/spline; 5=bell-spline)"
      write(lid,"(i15, a)") npar, " = npar:   Parameters number "
      write(lid,"(1pe15.8, a)") incl, &
         " = incl:   Inclination of wall in the convergent section [degrees]"
      write(lid,"(1pe15.8, a)") rcin, &
         " = rcin:   Radius of curvature at the entrance [m]"
      write(lid,"(1pe15.8, a, a)") rcgc, " = rcgc:   ", &
         "Radius of curvature in the throat in the convergent side [m]"
      write(lid,"(1pe15.8, a, a)") rcgd, " = rcgd:   ", &
         "Radius of curvature in the throat in the divergent side [m]"
      write(lid,"(1pe15.8, a)") rg, " = rg:     Throat radius [m]"
      write(lid,"(1pe15.8, a)") lch, " = lch:    Length of the chamber [m]"
      write(lid,"(1pe15.8, a)") lc, &
         " = lc:     Length of the convergent section [m]"

      do i = 1, npar

         if (ns == 2 .and. i == 2) then

            write(lid, "(1pe15.8, a, i1, a)") par(2)/pi*180.d0, &
               " = par(", i, "):"

         else

            write(lid, "(1pe15.8, a, i1, a)") par(i), " = par(", i, "):"

         end if

      end do

      write(lid,"(1pe15.8, a)") area, &
         " = area:   Area da parede da seccao divergente da tubeira [m2]"
      write(lid,*)

      if (area .ne. area) then

         write(*,*)
         write(*,*) "ERROR: error in calculus of the area."
         write(*,*)

         if (optm /= 0) call depp_save_fitness(0.d0, 5, &
            "ERROR: error in calculus of the area")

         stop

      end if

      prt = 1

      if (ns == 3) y(1) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, &
         rcin, x1, incl, x(1))

      if (ns == 4 .or. ns == 5) y(1) = bell_contour(lid, prt, rg, rcgc, rcgd, &
         rcin, x1, x5, incl, spl, ncp+2, xj, x(1))

      flush(lid)

      ! Write boundary coordinates

      open(10, file = trim(folder_output) // trim(sim_id) // "-boundary" // &
         nmbr // ".dat")

      write(10,*) "south boundary"
      do i = 1, nx-1
         write(10,*) x(i)-x0, 0.d0
      end do

      write(10,*) "north boundary"
      do i = 1, nx-1
         write(10,*) x(i)-x0, y(i)
      end do

      close(10)

   end subroutine get_boundary_nodes

   !============================================================================

   subroutine get_length_of_the_nozzle(lid, prt, optm, ns, npar, ncp, rg, &
         rcgc, rcgd, rcin, x0, x1, x5, incl, area, xj, spl, par)
      implicit none
      integer, intent(in) :: lid    !< Listing file id
      integer, intent(in) :: prt    !< Print nozzle contour coefficients (0=no; 1=yes)
      integer, intent(in) :: optm   !< Coupling with optimization method (0=no; 1=yes)
      integer, intent(in) :: ns     !< Nozzle shape (1=arbitrary; 2=conical; 3=parabolic; 4=bell-arc/spline; 5=bell-spline)
      integer, intent(in) :: npar   !< Number of parameters
      integer, intent(in) :: ncp    !< Number of control points
      real(8), intent(in) :: rg     !< Radius of the throat
      real(8), intent(in) :: rcgc   !< Radius of curvature in the throat in the convergent side [m]
      real(8), intent(in) :: rcgd   !< Radius of curvature in the throat in the divergent side [m]
      real(8), intent(in) :: rcin   !< Curvature radius in the input
      real(8), intent(in) :: x0     !< Coordinated at the entrance of the chamber [m]
      real(8), intent(in) :: x1     !< X coordinate of the entrance of convergent section
      real(8), intent(in) :: x5     !< X coordinate of the intersection of arc of throat with the divergent wall
      real(8), intent(in) :: incl   !< Inclination of the straight wall of the convergent section
      real(8), intent(in) :: area   !< Area of the surface wall
      real(8), intent(in) :: xj(npar-2)      !< x coordinates of the control points
      real(8), intent(in) :: spl(npar-3,4)   !< Coefficients of the spline
      real(8), intent(inout) :: par(npar)    !< Parameters of the divergent section

      ! Auxliary variables
      integer :: i, m
      real(8) :: aux, xaux1, xaux2, yaux1, yaux2, h, dist, raio, mult
      real(8), parameter :: pi = acos(-1.d0)

      m = 10000

      h = par(1)/m

      xaux1 = 0.d0

      select case (ns)

         case (1)

            yaux1 = arbitrary_contour(x0, x1, xaux1, incl, rcin, rcgc, rcgd, &
               rg, par)

         case (2)

            yaux1 = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, xaux1)

         case (3)

            yaux1 = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, x1, &
               incl, xaux1)

         case (4,5,6)

            yaux1 = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
               spl, ncp+2, xj, xaux1)

         case default

            write(*,*) "ERROR(1): invalid value for NS."

            stop

      end select

      aux = 0.d0

      do i = 1, 10*m

         xaux2 = i*h

         select case (ns)

            case (1)

               yaux2 = arbitrary_contour(x0, x1, xaux2, incl, rcin, rcgc, &
                  rcgd, rg, par)

            case (2)

               yaux2 = conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, &
                  xaux2)

            case (3)

               yaux2 = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, &
                  x1, incl, xaux2)

            case (4,5,6)

               yaux2 = bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, &
                  incl, spl, ncp+2, xj, xaux2)

            case default

               write(*,*) "ERROR(1): invalid value for NS."

               stop

         end select

         dist = dsqrt((xaux1 - xaux2)**2 + (yaux1 - yaux2)**2)

         raio = 0.5d0*(yaux1 + yaux2)

         aux = aux + 2.d0*pi*raio*dist

         if (aux >= area) then

            mult = (aux - area)/(2.d0*pi*raio*dist)

            xaux2 = xaux2 - mult*h

            if (ns == 3) then

               par(3) = parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, &
                  x1, incl, xaux2)

            end if

            par(1) = xaux2

            exit

         end if

         xaux1 = xaux2

         yaux1 = yaux2

      end do

      if (i == 10*m+1) then

         write(*,*)
         write(*,*) " ======== ERROR: It was not possible to obtain ========"
         write(*,*) " ======== the surface area of ​​wall required.   ========"
         write(*,*)

         if (optm /= 0) then

            call depp_save_fitness(0.d0, 4, "ERROR: not obtained surface area")

         end if

         stop

      end if

   end subroutine get_length_of_the_nozzle

   !============================================================================

   real(8) function arbitrary_contour(x0, x1, xi, incl, rcin, rcgc, rcgd, rg, &
         par)
      implicit none
      real(8), intent(in) :: x0     !< x coordinate of the chamber entrance
      real(8), intent(in) :: x1     !< x coordinate of the entrance of the convergente section
      real(8), intent(in) :: xi     !< x coordinate of interest
      real(8), intent(in) :: incl   !< Inclination of the straight wall of the convergent section
      real(8), intent(in) :: rcin   !< Curvature radius in input
      real(8), intent(in) :: rcgc   !< Radius of curvature in the throat in the convergent side [m]
      real(8), intent(in) :: rcgd   !< Radius of curvature in the throat in the divergent side [m]
      real(8), intent(in) :: rg     !< Radius of the throat
      real(8), intent(in) :: par(1) !< Lenght of de divergent section

      real(8) :: yi     !< y Coordinate of interest
      real(8) :: aux
      real(8), parameter :: pi = acos(-1.d0)

      aux = x0 + x1 + xi + incl + rcin + rcgc + rcgd + rg + par(1) ! calculation to avoid warnings

      yi = incl*xi**2 + rg

      arbitrary_contour = yi

   end function arbitrary_contour

   !============================================================================

   real(8) function conical_contour(par, rg, rcgc, rcgd, rcin, x1, incl, xi)
      implicit none
      real(8), intent(in) :: rg     !< Radius of the throat
      real(8), intent(in) :: rcgc   !< Radius of curvature in the throat in the convergent side [m]
      real(8), intent(in) :: rcgd   !< Radius of curvature in the throat in the divergent side [m]
      real(8), intent(in) :: rcin   !< Curvature radius in the input
      real(8), intent(in) :: x1     !< x coordinate of the entrance of convergent section
      real(8), intent(in) :: incl   !< Inclination of the straight wall of the convergent section
      real(8), intent(in) :: xi     !< x coordinate of interest
      real(8), intent(in) :: par(2) !< Parameters of the divergent section

      real(8) :: yi     !< y coordinate of interest
      real(8) :: x2     !< x coordinate of the intersection of 1st arc with the convergent wall
      real(8) :: y2     !< y coordinate of the intersection of 1st arc with the convergent wall
      real(8) :: x3     !< x coordinate of the intersection of convergent wall with the arc of throat
      real(8) :: y3     !< y coordinate of the intersection of convergent wall with the arc of throat
      real(8) :: x5     !< x coordinate of the intersection of arc of throat with the divergent wall
      real(8) :: y5     !< y coordinate of the intersection of arc of throat with the divergent wall
      real(8) :: xc1    !< x coordinate of the center of circumference of the curvature of entrance
      real(8) :: yc1    !< y coordinate of the center of circumference of the curvature of entrance
      real(8) :: xc4    !< x coordinate of the center of circumference of the curvature of throat
      real(8) :: yc4c   !< y coordinate of the center of circumference of the curvature of throat in the convergent side
      real(8) :: yc4d   !< y coordinate of the center of circumference of the curvature of throat in the divergent side
      real(8) :: ac     !< coefficient of the straight wall of the convergente section
      real(8) :: bc     !< Coefficient of the straight wall of the convergente section
      real(8) :: ad     !< Coefficient of the straight wall of the divergent section
      real(8) :: bd     !< Coefficient of the straight wall of the divergent section
      real(8), parameter :: pi = acos(-1.d0)

      xc4 = 0.d0
      yc4c = rcgc + rg
      yc4d = rcgd + rg

      ! straight wall of the convergent section
      bc = -tan(incl*pi/180.d0)
      x3 = -sqrt(bc**2*rcgc**2/(1 + bc**2))
      y3 = -sqrt(rcgc**2 - x3**2) + yc4c
      ac = -x3*bc + y3

      ! circle arc between the chamber and the convergent section
      x2 = x1 + sqrt(bc**2*rcin**2/(1 + bc**2))
      y2 = ac + bc*x2
      yc1 = y2 - sin((90.d0 - incl)*pi/180.d0)*rcin
      xc1 = x1

      ! straight wall of the divergent section
      bd = tan(par(2))
      x5 = sqrt(bd**2*rcgd**2/(1 + bd**2))
      y5 = -sqrt(rcgd**2 - x5**2) + yc4d
      ad = -x5*bd + y5

      if (xi < x1) then

         yi = yc1 + rcin

      else if (xi < x2) then

         yi = sqrt(rcin**2 - (xi - xc1)**2) + yc1

      else if (xi < x3) then

         yi = ac + bc*xi

      else if (xi < xc4) then

         yi = -sqrt(rcgc**2 - xi**2) + yc4c

      else if (xi < x5) then

         yi = -sqrt(rcgd**2 - xi**2) + yc4d

      else
         yi = ad + bd*xi

      end if

      conical_contour = yi

   end function conical_contour

   !============================================================================

   real(8) function parabolic_contour(lid, prt, par, rg, rcgc, rcgd, rcin, x1, &
         incl, xi)
      implicit none
      integer, intent(in) :: lid    !< Listing file id
      integer, intent(in) :: prt    !< Print nozzle contour coefficients (0=no; 1=yes)
      real(8), intent(in) :: rg     !< Radius of the throat
      real(8), intent(in) :: rcgc   !< Radius of curvature in the throat in the convergent side [m]
      real(8), intent(in) :: rcgd   !< Radius of curvature in the throat in the divergent side [m]
      real(8), intent(in) :: rcin   !< Curvature radius in the input
      real(8), intent(in) :: x1     !< x coordinate of the entrance of convergent section
      real(8), intent(in) :: incl   !< Inclination of the straight wall of the convergent section
      real(8), intent(in) :: par(3) !< Parameters of the divergent section
      real(8), intent(in) :: xi     !< x coordinate of interest

      real(8) :: yi     !< y coordinate of interest
      real(8) :: x2     !< x coordinate of the intersection of 1st arc with the convergent wall
      real(8) :: y2     !< y coordinate of the intersection of 1st arc with the convergent wall
      real(8) :: x3     !< x coordinate of the intersection of convergent wall with the arc of throat
      real(8) :: y3     !< y coordinate of the intersection of convergent wall with the arc of throat
      real(8) :: x5     !< x coordinate of the intersection of arc of throat with the divergent wall
      real(8) :: y5     !< y coordinate of the intersection of arc of throat with the divergent wall
      real(8) :: x6     !< x coordinate of the nozzle exit
      real(8) :: y6     !< y coordinate of the nozzle exit
      real(8) :: dy5    !< Inclination of wall in the intersection between the arc and the parable
      real(8) :: xc1    !< x coordinate of the center of circumference of the curvature of entrance
      real(8) :: yc1    !< y coordinate of the center of circumference of the curvature of entrance
      real(8) :: xc4    !< x coordinate of the center of circumference of the curvature of throat
      real(8) :: yc4c   !< y coordinate of the center of circumference of the curvature of throat in the convergent side
      real(8) :: yc4d   !< y coordinate of the center of circumference of the curvature of throat in the divergent side
      real(8) :: ac     !< Coefficient of the straight wall of the convergente section
      real(8) :: bc     !< Coefficient of the straight wall of the convergente section
      real(8) :: a      !< Coefficient of the parabolic wall of the divergent section
      real(8) :: b      !< Coefficient of the parabolic wall of the divergent section
      real(8) :: c      !< Coefficient of the parabolic wall of the divergent section
      real(8) :: det0
      real(8), parameter :: pi = acos(-1.d0)

      xc4 = 0.d0
      yc4c = rcgc + rg
      yc4d = rcgd + rg

      ! straight wall of the convergent section
      bc = -tan(incl*pi/180.d0)
      x3 = -sqrt(bc**2*rcgc**2/(1 + bc**2))
      y3 = -sqrt(rcgc**2 - x3**2) + yc4c
      ac = -x3*bc + y3

      ! circle arc between the chamber and the convergent section
      x2 = x1 + sqrt(bc**2*rcin**2/(1 + bc**2))
      y2 = ac + bc*x2
      yc1 = y2 - sin((90.d0 - incl)*pi/180.d0)*rcin
      xc1 = x1

      ! parabolic wall of the divergent section
      x5 = par(2)
      y5 = -sqrt(rcgd**2 - (x5 - xc4)**2) + yc4d
      dy5 = (x5 - xc4)/sqrt(rcgd**2 - (x5 - xc4)**2)
      x6 = par(1)
      y6 = par(3)

      det0 = 2*dy5*y5**2 + dy5*y6**2 - 2*dy5*y5*y6 - dy5*y5**2
      a = (y5 + dy5*x6 - y6 - x5*dy5)/det0
      b = (2*dy5*y5*x5 + y6**2 - 2*dy5*y5*x6 - y5**2)/det0
      c = (2*dy5*y5**2*x6 + dy5*y6**2*x5 + y5**2*y6 - y5*y6**2 &
         - 2*dy5*y5*y6*x5 - dy5*y5**2*x6)/det0

      if (xi < x1) then

         yi = yc1 + rcin

      else if (xi < x2) then

         yi = sqrt(rcin**2 - (xi - xc1)**2) + yc1

      else if (xi < x3) then

         yi = ac + bc*xi

      else if (xi < xc4) then

         yi = -sqrt(rcgc**2 - xi**2) + yc4c

      else if (xi < x5) then

         yi = -sqrt(rcgd**2 - xi**2) + yc4d

      else

         yi = (-b + sqrt(b**2 - 4*a*(c - xi)))/(2*a)

      end if


      if (prt == 1) then

         write(lid,"(1pe15.8, a)") a, &
            " = coefa:  Coefficient of nozzle parabolic contour"

         write(lid,"(1pe15.8, a)") b, &
            " = coefb:  Coefficient of nozzle parabolic contour"

         write(lid,"(1pe15.8, a)") c, &
            " = coefc:  Coefficient of nozzle parabolic contour"

      end if

      parabolic_contour = yi

   end function parabolic_contour

   !============================================================================

   real(8) function bell_contour(lid, prt, rg, rcgc, rcgd, rcin, x1, x5, incl, &
         spl, npar, xj, xi)
      implicit none
      integer, intent(in) :: lid    !< Listing file id
      integer, intent(in) :: prt    !< Print nozzle contour coefficients (0=no; 1=yes)
      integer, intent(in) :: npar   !< Number of parameters
      real(8), intent(in) :: rg     !< Radius of the throat
      real(8), intent(in) :: rcgc   !< Radius of curvature in the throat in the convergent side [m]
      real(8), intent(in) :: rcgd   !< Radius of curvature in the throat in the divergent side [m]
      real(8), intent(in) :: rcin   !< Radius of the curvature in the input
      real(8), intent(in) :: x1     !< x coordinate of the entrance of convergent section
      real(8), intent(in) :: x5     !< x coordinate of the intersection of arc of throat with the divergent wall
      real(8), intent(in) :: incl   !< Inclination of the straight wall of the convergent section
      real(8), intent(in) :: xi     !< x coordinate of interest
      real(8), intent(in) :: xj(npar-2)      !< x coordinates of the control points
      real(8), intent(in) :: spl(npar-3,4)   !< Coefficients of the spline

      integer :: i
      integer :: j
      integer :: k
      real(8) :: yi     !< y coordinate of interest
      real(8) :: x2     !< x coordinate of the intersection of 1st arc with the convergent wall
      real(8) :: y2     !< y coordinate of the intersection of 1st arc with the convergent wall
      real(8) :: x3     !< x coordinate of the intersection of convergent wall with the arc of throat
      real(8) :: y3     !< y coordinate of the intersection of convergent wall with the arc of throat
      real(8) :: xc1    !< x coordinate of the center of circumference of the curvature of entrance
      real(8) :: yc1    !< y coordinate of the center of circumference of the curvature of entrance
      real(8) :: xc4    !< x coordinate of the center of circumference of the curvature of throat
      real(8) :: yc4c   !< y coordinate of the center of circumference of the curvature of throat in the convergent side
      real(8) :: yc4d   !< y coordinate of the center of circumference of the curvature of throat in the divergent side
      real(8) :: ac     !< Coefficient of the straight wall of the convergente section
      real(8) :: bc     !< Coefficient of the straight wall of the convergente section
      real(8), parameter :: pi = acos(-1.d0)

      xc4 = 0.d0
      yc4c = rcgc + rg
      yc4d = rcgd + rg

      ! straight wall of the convergent section
      bc = -tan(incl*pi/180.d0)
      x3 = -sqrt(bc**2*rcgc**2/(1.d0 + bc**2))
      y3 = -sqrt(rcgc**2 - x3**2) + yc4c
      ac = -x3*bc + y3

      ! circle arc between the chamber and the convergent section
      x2 = x1 + sqrt(bc**2*rcin**2/(1.d0 + bc**2))
      y2 = ac + bc*x2
      yc1 = y2 - sin((90.d0 - incl)*pi/180.d0)*rcin
      xc1 = x1

      if (xi < x1) then       ! nozzle chamber

         yi = yc1 + rcin

      else if (xi < x2) then  ! arc of the entrance

         yi = sqrt(rcin**2 - (xi - xc1)**2) + yc1

      else if (xi < x3) then  ! straight wall of the convergent section

         yi = ac + bc*xi

      else if (xi < xc4) then  ! arc of the throat

         yi = -sqrt(rcgc**2 - xi**2) + yc4c

      else if (xi < x5) then  ! arc of the throat

         yi = -sqrt(rcgd**2 - xi**2) + yc4d

      else                    ! wall of the divergent section

         do k = 1, npar-4

            if (xj(k+1) >= xi) exit

         end do

         yi = spl(k,1)*(xi - xj(k+1))**3 + spl(k,2)*(xi - xj(k+1))**2 &
            + spl(k,3)*(xi - xj(k+1)) + spl(k,4)

      end if

      if (prt == 1) then

         do i = 1, npar-3

            write(lid,"(1pe15.8, a, i1, a)") xj(i), " = xj(", i, &
               "):  Control point of nozzle shape"

            do j = 1, 4

               write(lid,"(1pe15.8, a, i1, a, i1, a)") spl(i,j), " = spl(", i, &
                  ",", j, "): Spline coefficients of nozzle shape"

            end do

         end do

         write(lid,"(1pe15.8, a, i1, a)") xj(i), " = xj(", i, &
            "):  Control point of nozzle shape"

      end if

      bell_contour = yi

   end function bell_contour

   !============================================================================

   subroutine get_spline_coefficients(y5, dy5, Iexit, xj, yaux, n, spl)
      implicit none
      integer, intent(in) :: n         !< Number of equations of the spline
      real(8), intent(in) :: y5        !< y coordinate of the intersection of arc of throat with the divergent wall
      real(8), intent(in) :: dy5       !< Inclination of wall in the intersection between the arc and the spline
      real(8), intent(in) :: Iexit     !< Inclination of wall in the nozzle exit [degrees]
      real(8), intent(in) :: xj(0:n)   !< x coordinates of the control points
      real(8), intent(in) :: yaux(n)   !< y coordinates of the control points
      real(8), intent(out) :: spl(n,4) !< Coefficients of the spline

      integer :: i
      real(8), parameter :: pi = dacos(-1.d0)
      real(8) :: dye
      real(8) :: yj(0:n)
      real(8) :: b(0:n)
      real(8) :: g(0:n)
      real(8) :: h(n)
      real(8) :: a(0:n,0:n)

      dye = Iexit*pi/180.d0

      dye = dsin(dye)/dcos(dye)

      yj(0) = y5
      do i = 1, n
         yj(i) = yaux(i)
         h(i) = xj(i) - xj(i-1)
      end do

      a = 0.d0

      a(0,0) = -h(1)*0.3333333333333333d0
      a(0,1) = -h(1)*0.1666666666666666d0
      b(0) = (yj(0) - yj(1))/h(1) + dy5

      if (Iexit > 180.d0) then ! linear extrapolation in the throat
         a(0,0) = 1.d0
         a(0,1) = 0.d0
         b(0) = 0.d0
      end if

      do i = 1, n-1
         a(i,i) = 2.d0*(h(i) + h(i+1))
         a(i,i-1) = h(i)
         a(i,i+1) = h(i+1)
         b(i) = 6.d0*((yj(i+1) - yj(i))/h(i+1) - (yj(i) - yj(i-1))/h(i))
      end do

      if (Iexit > 90.d0) then
         a(n, n) = 1.d0
         b(n) = 0.d0
      else
         a(n,n-1) = h(n)*0.1666666666666666d0
         a(n,n) = h(n)*0.3333333333333333d0
         b(n) = (yj(n-1) - yj(n))/h(n) + dye
      end if

      g = 0.d0
      call tdma(a, b, g, n+1)

      do i = 1, n
         spl(i,1) = (g(i) - g(i-1))/(6*h(i))
         spl(i,2) = g(i)/2
         spl(i,3) = (yj(i) - yj(i-1))/h(i) + (2*h(i)*g(i) + g(i-1)*h(i))/6
         spl(i,4) = yj(i)
      end do

   end subroutine get_spline_coefficients

   !============================================================================

   ! Rescales the coefficients a and source of a linear system of 5 diagonals
   ! in order to avoid rouding errors.
   subroutine get_a5d_b_rescaling(nx, ny, a, b)
      implicit none
      integer, intent(in) :: nx     !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny     !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(inout) :: a(nx*ny,5)   !< Coefficients of the linear system
      real(8), intent(inout) :: b(nx*ny)     !< Source of the linear system

      integer :: i, j, np, npw, npe, nps, npn, npsw, npse, npnw, npne, iaux
      real(8) :: F

      ! West boundary
      i = 1
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npe = np + 1
         F = a(npe,3)/a(np,3)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         F = a(npw,3)/a(np,3)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! South boundary
      j = 1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         npn = np + nx
         F = a(npn,3)/a(np,3)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! North boundary
      j = ny
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         nps = np - nx
         F = a(nps,3)/a(np,3)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! Corners

      ! SW
      i = 1
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npne = npn + 1
      F = a(npne,3)/a(np,3)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

      ! SE
      i = nx
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npnw = npn - 1
      F = a(npnw,3)/a(np,3)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

      ! NW
      i = 1
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npse = nps + 1
      F = a(npse,3)/a(np,3)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

      ! NE
      i = nx
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npsw = nps - 1
      F = a(npsw,3)/a(np,3)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

   end subroutine get_a5d_b_rescaling

   !============================================================================

   ! Rescales the coefficients a and source of a linear system of 9 diagonals
   ! in order to avoid rouding errors.
   subroutine get_a9d_b_rescaling(nx, ny, a, b)
      implicit none
      integer, intent(in) :: nx  !< Number of volumes in csi direction (real + fictitious)
      integer, intent(in) :: ny  !< Number of volumes in eta direction (real + fictitious)
      real(8), intent(inout) :: a(nx*ny,9)   !< Coefficients of the linear system
      real(8), intent(inout) :: b(nx*ny)     !< Source of the linear system

      integer :: i, j, np, npw, npe, nps, npn, npsw, npse, npnw, npne, iaux
      real(8) :: F

      ! West boundary
      i = 1
      do j = 2, ny-1
         np   = nx*(j - 1) + i
         npe  = np + 1
         F = a(npe,5)/a(np,5)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! East boundary
      i = nx
      do j = 2, ny-1
         np = nx*(j - 1) + i
         npw = np - 1
         F = a(npw,5)/a(np,5)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! South boundary
      j = 1
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np = iaux + i
         npn = np + nx
         F = a(npn,5)/a(np,5)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! North boundary
      j = ny
      iaux = nx*(j - 1)
      do i = 2, nx-1
         np   = iaux + i
         nps  = np - nx
         F = a(nps,5)/a(np,5)
         b(np) = b(np)*F
         a(np,:) = a(np,:)*F
      end do

      ! Corners

      ! SW
      i = 1
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npne = npn + 1
      F = a(npne,5)/a(np,5)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

      ! SE
      i = nx
      j = 1
      np = nx*(j - 1) + i
      npn = np + nx
      npnw = npn - 1
      F = a(npnw,5)/a(np,5)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

      ! NW

      i = 1
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npse = nps + 1
      F = a(npse,5)/a(np,5)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

      ! NE
      i = nx
      j = ny
      np = nx*(j - 1) + i
      nps = np - nx
      npsw = nps - 1
      F = a(npsw,5)/a(np,5)
      b(np) = b(np)*F
      a(np,:) = a(np,:)*F

   end subroutine get_a9d_b_rescaling

   !============================================================================

   subroutine convert_int_to_char3(int3, char3)
      implicit none
      character*3 char3
      integer int3

      if (int3 < 10) then
         write(char3,"('00',i1)") int3
      else if (int3 < 100) then
         write(char3,"('0',i2)") int3
      else
         write(char3,"(i3)") int3
      end if

   end subroutine convert_int_to_char3

   !============================================================================

end module user
